Thermally-driven flows and turbulence in vibrated liquids

Within the vast array of variants encompassed by thermal (buoyancy) convection, the particular interest of the present work is on the specific dynamics which are enabled when standard steady gravity is replaced by a time-periodic body force induced by vibrations. The study is designed as a set of separate problems, where each exemplar aims to unravel the implications of the fundamental properties of this type of flow. These include the symmetry of the emerging pattern as perceived by a real observer and as seen in a “ time-averaged space ” , the synchronous or non-synchronous response of the velocity field to the applied forcing, the magnitude of the so-called Thermofluid-dynamic (TFD) distortions and the peculiar route of evolution towards chaos. A kaleido- scope of previously unknown solutions is reported giving emphasis to some still poorly known aspects such as the complex nature of the textural transitions that take place in the flow as the Gershuni number is increased (from 3.30 × 10 2 to 5.00 × 10 7 for Pr = 15). It is shown that the low-frequency regime is relatively stable over this range. In addition to the standard quadrupolar pattern, in such a case peculiar convective structures emerge where the time-averaged rolls display a very regular columnar arrangement, which has been rarely observed in earlier studies. Chaotic states are enabled when larger frequencies of vibration are considered. While for in- termediate frequencies concurrent aspects of the Feigenbaum and Manneville and Pomeau mechanisms can be recognized, the hallmark of the high frequency regime is its adherence to the standard Ruelle-Takens scenario.


Introduction
Since its beginnings, convection of natural origin hasn't ceased to amaze us with its peculiarity and counter-intuitive behaviors. Since the seminal experiments conducted by Lord Rayleigh more than one hundred years ago [1,2], fascination with the spatiotemporal patterns that arise spontaneously when non-isothermal fluids are subjected to forces of various natures has stimulated systematic experimental and theoretical research. All these efforts have finally resulted in the establishment of a specific area of study, nowadays generally known as the field of 'thermal convection'. Countless related technological applications can be found in many fields. Moreover, this form of convection is highly relevant to various natural sciences. From a purely theoretical standpoint, the main objectives of the related lines of research have been the analysis of the patterns of fluid motion in the course of evolution from a given initial state as some relevant (characteristic) parameters are changed. Stripped down to their fundamentals, these problems have proved themselves to be a source of complex non-linear effects leading to a myriad of possible manifestations and variants. These, in turn, have attracted the interest of investigators with various backgrounds and perspectives, encouraging them to elaborate a peculiar way of thinking and specific strategies of attack.
As an example of this modus operandi, it is worth recalling that the specific approach adopted so far by the community of theoretical physicists and engineers, is essentially based on the paradigm of a series of different canonical geometries. This practice has been motivated by the need to study the evolution of fluid flow under specific symmetries and discern how these symmetries are inherited or violated by the flow as the aforementioned characteristic numbers are varied. This specific methodology has its theoretical foundation in the intention to extract from specific cases, general features that can be used in the description of more complex configurations. This explains why for each geometry, attempts have been made to build a precise map in the space of parameters linking the possible input conditions to the 'output'.
Along these lines, it has been well understood that these problems consistently bare instabilities, which require dedicated analysis with the support of adequate mathematical tools to reveal the underlying 'physical' mechanisms. In addition to the "shape" of the cavity, crucial to these is the specific force driving the fluid flow.
As a fleeting glimpse into the literature immediately reveals, many studies have been based around the very classical two-dimensional (2D) square cavity. This configuration has become over the years a prototypical benchmark for numerical simulations. This is due to its simplicity and the high number of intrinsic symmetries (reflectional symmetries with respect to the vertical, horizontal and diagonal directions) that can be violated even though the tendency towards three-dimensionality is neglected. For the case of steady gravity and horizontal temperature difference (differentially heated cavity with vertical isothermal walls), many interesting numerical analyses have been produced for different values of the Prandtl number (ratio of the fluid kinematic viscosity and thermal diffusivity, Pr=ν/α) [3][4][5][6][7][8]; it is known that the flow can even develop turbulence [9][10][11][12][13][14][15][16][17][18][19][20]. These findings have been instrumental in understanding the role played by the (square) shape of the cavity and revealing the possible modes of convection and their peculiar path of progression as the control parameter (the Rayleigh number) is increased (steady→periodic→quasi-periodic→turbulent flow).
Despite all these valuable contributions, our knowledge of these phenomena is still relatively limited if one considers the equivalent circumstance in which fluid flow is produced by a time periodic body force in place of a steady body force.
A heuristic way to justify this substantial lack of knowledge is to consider the scarce availability of orbiting platforms (compared to terrestrial conditions), which has hindered for a relatively long time the possibility to conduct experiments in weightless conditions. In recent years, this specific variant of buoyancy flow, known as "thermovibrational convection", has started to receive significant attention with the advent of microgravity science and the extensive exploitation of the International Space Station [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35]. Interestingly, many of these studies have been motivated by the need to evaluate a priori the effects induced by the so-called "g-jitters". These are the inertial disturbances produced by pulsating or impulsive external loads typical of the considered orbiting platform such as motor firings during maneuvers, crew activity, station motorized equipments, station moving elements, space debris impacts, and the station structural elastic response to these loads [36][37][38][39][40][41][42][43][44][45]. A large cross section of fundamental research on this topic has been reviewed in the books by Gershuni and Lyubimov [21] and Lappa [46][47][48].
By virtue of these investigations, only disjoint glimpses of qualitatively and quantitatively different results in different sub-regions of the parameter space have been made available. Nevertheless, some important features or general characteristics of thermovibrational convection have been recognized. For example, it has been understood that a clear distinction (a kind of dichotomy) must be invoked between regimes where the vibrations have large amplitude and small frequency and vice versa.
In the former case, a direct cause-and-effect relationship can be established between the forcing and the ensuing flow. From a purely practical standpoint, this means that a non isothermal fluid generally reacts by producing a velocity field that is synchronous with the applied vibrations. This is equivalent to stating that, in analogy with the behavior of the driving force, by considering the average value of the velocity in a given point of the physical domain one would obtain a value that is almost exactly equal to zero (the average value taken by the acceleration resulting from the vibrations is indeed exactly 0 as we will show mathematically in Section 2).
However, this condition is no longer a property of the convective states that are obtained when the opposite condition with large frequency and small amplitude vibrations is considered. In these cases (Gershuni regime), the flow can develop a remarkable (non-negligible) time-averaged velocity component standing out from the instantaneous flow, which is essentially time-periodic [49][50][51][52][53][54].
Superimposed on this complexity is the ability that these flows (in both limits of small and high frequencies) to develop instabilities when the relevant control parameter (the vibrational Rayleigh number) is increased. These instabilities can contribute to make the frequency spectrum progressively more complex until a completely turbulent state is attained.
The present analysis aims to improve our understanding of this form of convection and its possible manifestations by integrating the information reported in earlier investigations and focusing on previously overlooked aspects. More specifically, assuming the classical differentially heated square cavity as reference case, the following questions are tackled here: the existence and the dynamics of regimes of regular fluid motion characterized by different degrees of spatial symmetry, the reverberation of these symmetries and their violation (or rupture) on the underlying temperature field (leading to the so-called thermofluid-dynamic "distortions"), the influence of the frequency of vibrations on such symmetries and finally, the peculiar path taken by this system in its progression towards chaos (an aspect still almost completely unknown for what concerns thermovibrational convection). As anticipated in the introduction, for simplicity and in line with a vast portion of the existing literature, we focus on a square differentially heated box with characteristic size L (see Fig. 1).
Microgravity conditions are considered, i.e. steady gravity is absent; rather, the force driving the fluid flow is a time-periodic form of buoyancy produced by the combined effect of an imposed temperature difference and a sinusoidal displacement of the cavity in time. As it can easily be demonstrated according to simple mathematical arguments, such a displacement can be considered equivalent to a time-periodic acceleration. Indeed, representing the vibrations via the relationship: where b is the amplitude, ω=2πf is the angular frequency of the displacement and n is the unit vector along the direction of vibrations, the ensuing time-varying acceleration can be expressed mathematically by taking the second derivative of Eq. (1) with respect to time, i.e.
Notably, eq. (2) satisfies the condition: i.e. (as anticipated in the introduction) its time-averaged value over one period of oscillation 2π/ω is exactly zero. As the reader might have realized at this stage, the most important implication of Eq. (3) is the indication that it provides about time-averaged effects in the fluid flow, which (if present) must necessarily be ascribed to the non-linear nature of the governing equations (presented in Section 2.2).

The governing equations
The non-dimensional balance equations for mass, momentum and energy can be written as: where β T is the thermal expansion coefficient (used to make the relationship between fluid density and temperature linear as required by the standard Boussinesq approximation) and ΔT is the imposed temperature difference. The non-dimensional form of these equations follows naturally from the choice of suitable reference quantities, these being in the present case: L, L 2 /α, α/L, ρα 2 /L 2 , ΔT for the Cartesian coordinates (x,y), time (t), velocity (V), pressure (p) and temperature (T), respectively. The parameter (bω 2 β T ΔTL 3 )/να = Ra ω appearing in Eq. (5) is an analogue of the classical Rayleigh number used in earlier studies dealing with classical buoyancy convection (steady gravity case). However, the problem also strongly depends on the relative orientation of the imposed temperature gradient and of the direction of the periodic acceleration [24,47]. In particular, as shown in Fig. 1, in the present work, n is perpendicular to the gradient of temperature(which may be regarded as a distinguishing mark with respect to the companion study by Crewdson and Lappa [55]).

Projected equations and system intrinsic symmetries
For the conditions shown in Fig. 1, the projection of the momentum equation on the coordinate axes simply reads: where u and v are the velocity components along x and y, respectively, and ( V⋅∇V Eqs. (7) and (8) are particularly instructive if they are considered together with the related kinematic and thermal boundary conditions. Indeed, they can be used to elaborate some initial very useful considerations about the possible symmetries of the emerging flow.
Following the work of Lappa [30][31][32][33] and Lappa and Burel [34], here we assume no-slip conditions for all the solid walls, while from a thermal point of view the walls parallel and perpendicular to the vibrations are isothermal and adiabatic, respectively. Moreover, in terms of initial conditions at the instant (t=0), we consider a linear temperature distribution along the x axis, i.e. T=x, while the fluid is postulated to be at rest (i.e. V=0 everywhere).
Taking into account this specific combination of directions of temperature gradient, boundaries and driving force, using the spatial symmetry as a classification criterion, the set of possible convective modes in a square cavity can be split into four main categories as follows: • (sa): The symmetric-antisymmetric mode. The number of emerging convective cells along the direction of the temperature gradient (m x ) is odd while their number along the direction of the driving force (the vibrations in the present case) is even (m y ). The resulting pattern takes a symmetric configuration with respect to a line parallel to the x axis cutting the cavity in two identical halves. • (as): The antisymmetric-symmetric mode. The symmetry properties are inverted with respect to the (sa) case. The number of emerging convective cells along the direction of the temperature gradient (m x ) is even, while their number along the direction of the driving force (m y ) is odd. The resulting pattern takes a symmetric configuration with respect to a line parallel to the y axis cutting the cavity in two identical halves. • (ss): The symmetric-symmetric mode. The number of rolls is even regardless of the considered fundamental direction. The resulting pattern is symmetric with respect to both the two perpendicular lines (parallel to the x and y axes, respectively), passing through the geometrical center of the cavity. The rolls are arranged in 'columns' if the number of rolls along x is larger than the number of rolls along y, while they reduce to a pattern with the central symmetry if their number along these two directions is the same (m x =m y ). • (aa): The antisymmetric-antisymmetric mode. No specific symmetry can be recognized in this case. The number of rolls is odd along both the direction of the temperature gradient and the imposed vibrations (the simplest possible realization of such a condition being a single roll pervasive throughout the square cavity, i.e. m x =m y =1).
Though this classification was originally introduced by Mizushima and Mizushima and Adachi [56,57] for standard Rayleigh-Bénard convection, it displays a good degree of generality, which makes it applicable also to other contexts, including the case of thermovibrational flow [54].

Time averaged flow and the Gershuni number
Though Mizushima's arguments are still valid for thermovibrational flow, however, a proper distinction must be considered between the symmetry properties of the aforementioned so-called time-averaged flow, emerging as the consequence of purely non-linear effects, and the fluctuating part of the flow, representing the 'linear response' of the vibrated cavity to the applied forcing, namely: where P = 2π/Ω and In the present study, these quantities are determined 'a posteriori' after evaluating the fluid velocity V via direct numerical solution of the governing equations in their complete time-dependent and non-linear form (as illustrated in Section 3). This procedure may be regarded as a relevant tool to deepen the interpretations of these phenomena. Indeed, it offers valuable advantages with respect to alternate formulations elaborated in the literature where the user is prevented from getting any relevant information about the behavior of T' and V' (such as the socalled Gershuni formulation, by which the time-averaged quantities can be obtained immediately using a potential function) [49][50][51][52].
Although we do not rely on the Gershuni's simplified balance equations, as a concluding remark for this section, we report the definition of the so-called Gershuni number (Gs), see Eq. (11), known for its ability to influence the magnitude and symmetry properties of the timeaveraged flow [30][31][32][33][34].

The projection method
In the present section, a detailed description of the 'projection' (or 'primitive-variable') method used for the integration of Eqs. (4)-(6), is provided. From a purely mathematical point of view, these equations and the related boundary conditions should be regarded as a hybrid set of parabolic-elliptic-hyperbolic equations to be solved for the determination of the main problem unknowns, namely V, p and T. In particular, the so-called projection numerical techniques are heavily reliant on the coupling between the first two of these primitive variables, i.e. pressure and velocity (which explains why methods of such a kind are simply indicated as 'pressure-solvers' in most of existing commercial or opensource CFD packages).
From a practical standpoint, the numerical problem is formulated and solved in a multistage form. More precisely, the associated timemarching algorithm is structured around the main objective of systematically integrating an alternate set of equations (equivalent to the original set of mass and momentum balance equations from a mathematical point of view, but not from a numerical standpoint, Gresho [58]) in order to obtain a velocity field that satisfies at the same time the incompressibility constraint and the momentum equation. Obviously, this approach also relies on physical and mathematical principles involving well-posedness theory, which (as shown at the end of this section) are key to the physically consistent and computationally reliable method derivation.
The first alternate equation is obtained by dropping the pressure term in the momentum equation, i.e.
This truncated equation is no longer mathematically equivalent to the original balance of momentum (which explains why the velocity V* is generally referred to as 'provisional' velocity field). The quantity V*, however, still retains some physical meaning. This can be immediately verified by applying the curl operator to all of the terms of Eq. (12). Such an artifice would indeed result in the balance equation for the 'effective' flow vorticity; vorticity indeed is known to be independent ofthe gradient of pressure, whose curl is always mathematically equal to zero (owing to an intrinsic property of the nabla operator). In a second step this velocity is corrected to recover completely its physical relevance, that is, to make sure that it also satisfies the incompressibility constraint. This is formally achieved introducing a second alternate equation where the provisional velocity field is linearly combined with the previously neglected gradient of pressure, that is where C is a constant and p still plays the role of an unknown. A third (and final) alternate equation to determine the otherwise unknown pressure is finally obtained by forcing Eq. (13) into the continuity equation, which yields and effective elliptic equation for p: This extra equation formally closes the problem from a numerical point of view if C is set equal to Δt, where Δt is the time integration step. As mentioned earlier, this splitting can be seen as a realization of precise theorems. As an example, Eq. (13) should be regarded as a direct application of the Helmholtz decomposition theorem (with V and p playing the role of solenoidal part and scalar function underpinning V*, respectively [59]). Well-posedness is also guaranteed by the Ladyzhenskaya theorem [60], which states that a vector field can be considered fixed if its vorticity, divergence and component perpendicular to the domain boundary are assigned (in this regard, it is worth noticing that, as explained before, the divergence and curl of the velocity field determined using Eqs. (12)- (14) are mathematically identical to those potentially following the solution of the 'original' set of equations). Moreover, the velocity component perpendicular to the boundary is set to zero (accordingly, the pressure equation can be solved with homogeneous Neumann conditions [61][62][63]).
The above arguments may be regarded as a general description of the philosophy on which all projection methods are based. In particular, in the present work, the required computations have been carried out using the so-called PISO (Pressure Implicit Split Operator) technique, a specific variant of the general projection methods available in OpenFoam, based on the original implementation introduced by Issa [64]. This realization is based on a collocated grid approach, i.e. all the primitive variables are defined in the center of the computational cells; accordingly, the velocity requires a specific treatment on the faces of the cell (in order to improve the coupling of velocity and pressure), where it is interpolated using the special scheme by Rhie and Chow [65]. The energy equation for the determination of the temperature has been solved in a segregated manner. The convective and diffusive terms appearing in both the momentum and energy equations have been dealt with through upwind and central differencing schemes, respectively, while backward differencing in time has been implemented for the unsteady terms. As a concluding remark for this section, we wish to point out that the elliptic pressure equation (Eq. (14)) has been integrated resorting to a Generalized Geometric-Algebraic Multi-Grid (GAMG) strategy.

Validation
To demonstrate the reliability of the numerical method presented in Section 3.1 (and the related implementation available in OpenFOAM), preliminary simulations have been compared with well established benchmarks available in the literature. In particular, the patterns reported by Lappa [46] and Monti and Savino [54] have been selected given the similarity between their configuration and that examined in the present work. Indeed, in those studies, the steady gravity field was replaced by a time periodic body force and set perpendicular to the temperature gradient. The outcomes of such simulations are presented in Fig. 2 for Pr = 15 and Gs = 3.30 × 10 2 and 3.30 × 10 4 , where the pattern of the time-averaged velocity field is shown. The classical transition from the so-called quadrupolar flow field to the inversional symmetry pattern, as the Gershuni number is increased, can be clearly recognized there.
For the sake of additional validation, and given the scarcity of relevant results in the literature, the solution method has been verified further through 'ad hoc' comparisons between the results provided by OpenFoam and those yielded by an alternate computational platform (namely, the same software implemented by Lappa [30] and Lappa and Burel [34]). Although, the same class of pressure-velocity methods discussed before is at the root of both implementations, they differ in regard to temporal schemes (which are implicit and explicit for OpenFOAM and Lappa's algorithm, respectively). Another distinguishing mark is represented by the level of coupling for velocity and pressure (this being based on a staggered arrangement of primitive variables in the latter software).
The outcomes of the focused comparisons conducted for a series of values of the vibrational Rayleigh number (namely, Ra ω = 1.08 × 10 5 , Ra ω = 4.47 × 10 5 , Ra ω = 8.56 × 10 5 , and Ra ω = 1.08 × 10 6 , for Pr = 15 and Ω = 9.38 × 10 2 ) are quantitatively substantiated in Fig. 3. The reported frequencies have been obtained through a Fast Fourier Transform (the X, Y and Z axes representing the number of frequencies detected in the spectrum, the frequency value and the corresponding vibrational Rayleigh number, respectively).
As evident in this figure, the two codes provide an almost perfect agreement in terms of frequency spectrum (the only minor difference relates to the slightly higher number of frequencies obtained with the inhouse code in one of the cases, for which, however, the related amplitude has been found to be almost negligible).

Grid refinement study
Given that the specific nature of the considered problem may potentially lead to the development of instabilities and even transition to turbulence in some conditions, special attention has been paid to build relevant grids, i.e. meshes able to capture properly all the spatial scales involved in the emerging phenomena.
In the absence of relevant (observational) information to properly address this point (due to the limited number of earlier investigations on the subject), such a study has been based on general criteria available in the literature for problems that share common aspects with the present one.
As a first requirement to be met, the need to have a sufficient resolution inside the thermal boundary layers has been considered. Though a dichotomy is often drawn between buoyancy convection induced by steady gravity or a time-periodic acceleration, indeed, these phenomena display some affinities. For relatively high values of the Prandtl number, both are prone to develop regions of limited thickness located in proximity to the solid boundaries, where temperature displays strong gradients or steep profiles. As reported by several authors (see, e.g., the Order of Magnitude Analysis by Russo and Napolitano [66]) the thickness of this layer can be estimated as: As shown, e.g., by Shishkina et al. [67] 'a priori' knowledge of the thickness of the thermal boundary can be turned into precise numerical requirements by expressing the number of cells needed to capture  properly the dynamics of the thermal boundary layer as a function of Ra, i.e.
Superimposed on these constraints is the need to capture the smallest flow structure eventually present in the fluid flow when it enters the turbulent regime. The typical size of these features (in general manifesting themselves in the form of 'eddies') can be evaluated through the so-called Kolmogorov length scale. For cavities undergoing a body force perpendicular to the imposed temperature gradient this quantity can be estimated using the correlation originally elaborated by Paolucci [10] and Farhangnia et al. [13].
A summary of the theoretical grids defined using the above criteria under the assumption of uniform mesh are presented in Table 1, together with the corresponding thickness of the thermal boundary layers (as provided by effective numerical simulations based on such grids). It can be noted that, in the frame of this initial assessment of the required spatial resolution, Eq. (16) has been found to be the most restrictive mesh-resolution controlling factor over the considered range of values of the vibrational Rayleigh number.
In order to circumvent difficulties deriving from the prohibitively expensive computing cost resulting from the enhanced computational overhead required by uniform grids satisfying all these requirements, we have decided to use a non-uniform distribution of computational cells.
Along these lines, Fig. 4 illustrates the different mesh types adopted for the ensuing additional mesh refinement study. It can be noted that a bias is included towards the walls of the cavity. By doing so, the number of nodes in the boundary layer is increased instead of increasing the number of nodes across the whole domain, thereby reducing the required computational power.
Though the above discussion might be regarded as an exhaustive treatment of the problem in terms of space resolution requirements, it does not consider another important restrictive condition that plays a crucial role when flows induced by non-constant driving forces are considered, i.e. the needed temporal resolution. This specific aspect requires the consideration of two additional concurrent factors.
For themovibrational flow, both the instantaneous and time average velocities can be analyzed by means of Fast Fourier Transforms (FFT), whereby the velocities are presented as a function of frequency rather than a function of time. The appropriate sampling rate must be chosen in order to ensure all frequencies are taken into account. Traditionally, the minimum sampling rate of a continuous signal must obey the Nyquist-Shannon Sampling Theorem. The theorem states that the signal must be sampled at a frequency superior to twice the frequency of the signal f s (the Nyquist Frequency), giving f s /2<f. The time step of the simulation is defined as: ∆t = 2π/ωNp, where, Np is the number of samples per period and must satisfy the Nyquist Sampling Theorem whereby Np >2f . The second constraint to be aware of is the well-known effect of the Courant number C (which must be smaller than one to ensure numerical algorithm stability). For all the simulations presented in this study, the time integration step (satisfying both abovementioned requirements) has been selected in such a way to have a minimum of 400 steps for each cycle of vibrational modulation.
In order to confirm the viability of all these theoretical requirements, an additional mesh refinement study has been carried out on two cases, shown in Fig. 5, where the frequency spectrum has been chosen as the comparative output parameter. In case (a) (Ω = 9.38 × 10 3 , Ra ω = 1.08 × 10 6 , Gs = 1.00 × 10 5 ), a grid resolution of 80 by 80 cells satisfies the theoretical requirements. Accordingly, two peripheral simulations have been carried out, one with a coarser grid (60 by 60 cells) and one with a finer grid (100 by 100 cells). The results from the higher resolution grid are quasi-identical to those obtained with the predicted 80 by 80 grid, therefore this grid resolution has been used in the present study for all cases with Ra ω <1.08 × 10 6 . As for case (b) (Ω = 9.38 × 10 3 , Ra ω = 7.66 × 10 6 , Gs = 5.00 × 10 6 ), a grid resolution of 125 by 125 cells satisfies the theoretical requirements. For this case, the two additional simulations have been carried out on a 100 by 100 cell grid for the coarser case, and a 140 by 140 cell grid for the finer case. It can be noted that for any mesh exceeding 100 by 100 cells a bias is applied at the walls as indicated in Fig. 4. As this case is considered weakly turbulent (this is expanded on in Section 4.1), the exact frequencies and amplitudes can vary from simulation to simulation, however, it is evident here that the response of the system is very similar for all three grid resolutions. Therefore, a grid size of 125 by 125 is acceptable in this case. These two representative cases have shown that the theoretical requirements based on the complex set of criteria illustrated before can be used as a relevant basis to determine 'a priori' the needed resolution.

Results
In the present section, the dynamics revealed by the numerical simulations are discussed with the aid and support of flow patterns showing both V (time-averaged flow) and the total velocity V. These serve for a better representation of the structure of convection in a 'timeaveraged space' and as it would be seen by a real observer, respectively.
Several focused examples and situations are considered, many of which are of a prototypical nature. In line with the companion study by Crewdson and Lappa [55], where the alternate configuration with vibrations parallel to ΔT was treated, the Prandtl number is set to Pr = 15, while a wide range of orders of magnitude is covered in terms of Rayleigh number (Ra ω ). The parametric investigation also includes assessment of the role played by Ω (allowed to span an interval that includes both limits of high and small frequencies).
To elaborate such a treatment, we start from the circumstances considered for the validation study, i.e. the two cases with Gs = 3.30 × 10 2 and 3.30 × 10 4 , respectively, and Ω = 9.38 × 10 3 (corresponding to the aforementioned high-frequency Gershuni regime). In the light of the arguments introduced in Sect. 2.3, these two patterns of time-averaged convection can immediately be interpreted in terms of 'fundamental modes' or combination of them. In particular, the case with Gs = 3.30 × 10 2 corresponds to the aforementioned archetype of quadrupolar flow field, i.e. a classical example of symmetric-symmetric (ss) mode with equal (even) number of rolls along the x and y directions (a rare example of flow with the central symmetry, m x =m y =2). When the Gershuni number is increased by two orders of magnitude, the convective structure with central symmetry is taken over by a diagonal mode, i.e. a flow displaying reflectional symmetry with respect to one of the diagonals of Table 1 Theoretical and measured boundary layer properties.

The high frequency case: the Ruelle-Takens scenario
In the present sub-section, the Gershuni number (vibrational Rayleigh number) is varied parametrically from the first reference case Gs = 3.30 × 10 2 (corresponding to Ra ω = 6.22 × 10 4 ) to Gs = 5.00 × 10 7 (Ra ω = 2.42 × 10 7 ), thereby spanning 5 orders of magnitude. Although a total of 12 cases are simulated for a given frequency Ω = 9.38 × 10 3 , in the following, only 6 of these cases are discussed to illustrate the different possible 'solutions' (the reader being referred to Table 2 for a complete description of all simulations and the respective regimes).
For 3.30 × 10 2 <Gs<5.00 × 10 3 (or 6.22 × 10 4 <Ra ω <2.40 × 10 5 ) the streamlines of the time-averaged velocity field show the quadrupolar flow field already discussed in relation to the first validation case (the first fundamental mode displaying central symmetry, see Fig. 6, left panel). Starting from a single roll, the instantaneous flow field (Fig. 6, right panels) changes periodically in time, whereby four small rolls emerge at the outer corners of the cavity and then disappear leaving behind them again one single roll (this behavior occurring twice per period). Following up on this point, it is worth recalling that the single vortex pervasive throughout the entire square cavity, regularly changes its sense of circulation (from the counterclockwise to the clockwise orientation) according to the sign of the acceleration.
As the order of magnitude of the Gershuni number shifts from O(10 3 ) to O(10 4 ), the time-averaged flow field transitions from the quadrupolar field to a diagonal state (the second fundamental mode, Fig. 7, left panel), in agreement with the validation case. The instantaneous flow field also undergoes some modifications; indeed, it displays a more involved sequence of snapshots, where occasional manifestations can be seen of a diagonal convective structure (t = t 0 +0.1P) due to the superposition of (ss) and (aa) modes and a columnar arrangement of rolls (m y =1, m x =3) embedded in a larger scale circulation for t= t 0 +0.8P.
On increasing the Gershuni number to Gs = 10 5 (Ra ω = 1.08 × 10 6 ), substantial changes can be spotted in both the time-averaged and instantaneous fields (Fig. 8). In what concerns the time-averaged velocity, the dominant diagonal orientation of the main circulation is conserved, however the two inner roll seems to be 'squashed' by the growing of two smaller rolls in the outer corners of the cavity. As for the instantaneous field, a time periodic scenario is still observed. However, an increase in complexity in the patterning behavior is evident. The single roll pervasive throughout the cavity ((aa) mode with m y =m x =1), seen for lower values of Ra ω and Gs, does no longer represent the preferred transitional stage of evolution; it is replaced by transitional   Moving on to the case Gs = 10 6 (Ra ω = 3.42 × 10 6 ), although the topological configuration of the time-averaged flow field remains unchanged, some morphological alterations can be seen (Fig. 9, left panel). The two inner rolls embedded into the main circulation, which previously were perfectly aligned with the NW-SE diagonal direction, are now displaced towards the lateral (external walls) and, at the same time, towards the cavity midplane x = 1/2 (thereby creating the illusion of an inclined Z-shaped vortex). This is accompanied again by an increase in complexity in the instantaneous field. As in this case the velocity signal is no longer strictly synchronous with the forcing frequency (this will be expanded on in the next section), a selection of snapshots illustrating the various possible multi-cellular states is presented in Fig. 9, right panel.
The next interesting transition occurs when Gs is increased to 5.00 × 10 6 (and Ra ω to 7.66 × 10 6 ). In such conditions, notably, the timeaveraged field presents a slight fluctuation in time (it becomes time dependent, Fig. 10a). Surprisingly, the 4-roll configuration is recovered. Qualitative comparison with Fig. 6, however, reveals that the distribution of streamlines is seemingly distorted.
For what concerns the instantaneous field (as the reader will easily realize by inspecting the related sequence of snapshots in Fig. 10b), the temporal evolution of the flow can be seen as a continuous swap between a two-roll diagonal mode (oriented along the NW-SE or the NE-   The next figure of the sequence simply illustrates that as Gs is increased to Gs = 5.00 × 10 7 (Ra ω = 2.42 × 10 7 ), both the time-averaged field and the instantaneous field become relatively chaotic ( Fig. 11a and   b, respectively). Symmetry along the diagonal axis is (mostly) lost in both the instantaneous and time averaged field. This is accompanied by an increase in the number of eddies that appear in the flow, emerging from multiple locations in the cavity.
In the following, the earlier arguments developed for the topology of  the streamlines of the time-averaged and instantaneous flow field and the related temporal and morphological evolution are complemented by a frequency spectrum analysis based on the velocity signal measured by a numerical probe located in the geometrical center of the cavity. To extract these frequency spectra from the velocity signal, for each simulation a Discrete Fourier Transform has been implemented using a Fast Fourier Transform algorithm (FFT), whereby the time domain is transformed to the frequency domain. In this regard, Fig. 12 is instrumental in revealing that (as expected) an increase in Gs for a fixed Ω (equivalent to an increase in Ra ω ) makes the number of detected signal frequencies higher. In the following, for simplicity we use the following symbols to indicate integer multiples of the forcing frequency present in the spectrum: Ω 2 = 2Ω, and Ω 3 = 3Ω. As the reader will realize through inspection of Fig. 12, over the extended range 3.30 × 10 2 <Gs<5.00 × 10 5 (6.22 × 10 4 <Ra ω <2.42 × 10 6 ), the emerging flow can be considered synchronous with the imposed forcing (as all detected frequencies are harmonics of the original non-dimensional frequency of vibrations Ω).
However, from Gs = 10 6 , Ra ω = 3.42 × 10 6 and Gs = 5 × 10 6 , Ra ω = 7.66 × 10 6 the FFTs show a wide range of frequencies. Additional frequencies are seen to appear at the right and left-hand side of Ω and its harmonics (Fig. 13). Although these frequencies are non-synchronous with the main frequencies Ω, Ω 2 and Ω 3 , they follow an interesting pattern whereby: Ω 2R = 2Ω 1R and similarly, 2Ω L = 2Ω 1L , where Ω 1R and Ω 1L represent incommensurate frequencies that appear to the right and left of Ω respectively.
As a final look at the lower panels in Fig. 13 would indicate, when the Gershuni number is increased beyond Gs = 5 × 10 6 (and the Rayleigh number beyond Ra ω = 7.66 × 10 6 ) the behavior becomes more turbulent (as witnessed by the continuous distribution of frequencies in the spectrum, which quantitatively substantiate the complex patterning behavior reported in Fig. 11).
From the streamlines of the instantaneous and time-averaged velocity fields, combined with the frequency spectrum analysis, the overall behavior of the considered system can therefore be categorized as illustrated in Table 2. As the appearance of turbulence (a continuous frequency spectrum is preceded by the existence of 3 incommensurate frequencies and related higher-order harmonics), we also argue that in the range of high frequencies (Ω≅ 10 4 ), the transition to chaos essentially adheres to the Ruelle-Takens-Newhouse scenario [69,70].

Intermediate frequency and the Feigenbaum sequence
Having completed a sketch of the dynamics in the high-frequency regime, in the present section we turn to considering circumstances for which the frequency of the imposed vibrations is smaller (decreased by one order of magnitude, i.e. Ω = 9.38 × 10 2 ).
Following the same approach already undertaken for the higher frequency case, first, a description of the velocity field in terms of its various manifestations as a function of Gs at fixed Ω is provided.
As evident in Fig. 14, for Gs = 10 5 , Ra ω = 1.08 × 10 5 , the time averaged field is steady in time and takes the form of a large inner roll accompanied by two very small rolls at two of the outer corners. The instantaneous field shows an evolving patterning behavior displaying various realizations and superposition of convective modes with (aa) and (as) symmetry (columnar arrangements of two or three rolls embedded into a larger scale circulation).

G. Crewdson and M. Lappa
On increasing the Gershuni number to Gs = 5.83 × 10 6 (Ra ω = 8.27 × 10 5 ), the dynamics become weakly chaotic as witnessed by the unsteady time averaged flow undergoing random switches between distorted three-roll and four-roll patterns. As shown in Fig. 15, the increased complexity of the flow is also very evident in the instantaneous velocity field, where, though the dominance of columnar modes can still be identified, a variety of textural transitions take place.
An interesting change occurs for 7.08 × 10 6 <Gs<1.00 × 10 7 , (9.12 × 10 5 <Ra ω <1.08 × 10 6 , see, e.g., Fig. 16), where, most surprisingly, a back transition to a steady state is observed in terms of time-averaged  field. Such a simplification can also be spotted in the instantaneous velocity, where only a limited number of textural transitions can be detected.
Setting Gs to 5.00 × 10 7 , Ra ω = 2.42 × 10 6 , however, makes the flow relatively turbulent (see Fig. 17). While the time-averaged field fluctuates between a three-roll structure similar to that seen in previous cases and more complex three, four and five-roll patterns in a seemingly random manner, the instantaneous velocity gives rise to a rich set of patterns, which include many eddies as already seen previously for chaotic flow in the high-frequency regime.
Given the counter-intuitive evolution displayed by the system as a function of the Gershuni number (or vibrational Rayleigh number), Fig. 16. Gs = 1.00 × 10 7 , Ra ω = 1.08 × 10 6 , Ω = 9.38 × 10 2 : time-averaged flow (left) and instantaneous flow (right), where the instantaneous snapshots are taken at 0.1P intervals from t 0 = 0P to t= t 0+ P. unlike Sect. 4.1(where the analysis was limited to determining the frequency spectrum and showing its increasing complexity with Gs), here we resort to a more sophisticated strategy. This (multi-faceted) approach involves additional tools generally used for the assessment or 'quantification' of the level of chaos embedded in a non-linear system.
To do so, inspiration is taken from the work by Ueno et al. [71] . Although this referenced study was focused on a different type of convection (Marangoni flow in liquid bridges), the authors demonstrated that an alternate approach (with respect to standard two-dimensional Poincarré maps) based on the reconstruction of the pseudo-phase space from the time series of the surface temperature can be used to acquire useful insights into the chaotic behavior of a flow. Although in the present work, the velocity time series is used in place of the surface temperature time series, the underlying principles are essentially the same. The sequence of plots in Fig. 18 represents the outcome of such an analysis. This figure illustrates in an ordered fashion the instantaneous velocity field, the related Fourier transform and finally the reconstructed phase space for increasing values of Gs (Ra ω ).
We wish to recall that this fruitful three-fold representation method has also been employed by other authors such as Paolucci and Chenworth [9] and Guzman and Amon [72]. In the present study, the phase space reconstruction has been achieved by using the time delay embedding procedure, stemming from the mathematical framework proposed by Takens [73]. Where the observed time series (in terms of v, where v is the velocity component along y) lies on an 'unobservable' attractor, the time series is used to reconstruct this attractor by plotting the lagged coordinated vectors: where m is the embedding dimension and τ is the time delay (or lag). This method is highly informative of the specific behavior of the dynamic system as an increase in complexity in the velocity signal leads to an increase in complexity of both the frequency spectrum and the reconstructed phase space. As a first example, looking at the panel 18(a), it can be seen that the velocity signal is initially monochromatic and strictly synchronous with the forcing; the associated reconstructed phase space is composed of a single loop.
When looking at panel 17(b) however, although the signal is still regular, the presence of a second frequency equal to Ω/2 can be recognized; accordingly, a second loop appears in the phase space. This in an indication that for Gs = 1.70 × 10 6 , Ra ω = 4.47 × 10 5 , a perioddoubling bifurcation has taken place.
Careful analysis of Fig. 17 (the reader being also referred to the information summarized in Table 3) also reveals that while for Gs = 5.83 × 10 6 , Ra ω = 8.27 × 10 5 (panel 18(c)) the flow can be considered turbulent, for Gs = 6.25 × 10 6 , Ra ω = 8.56 × 10 5 (panel 18(d)), a twofrequency state (consisting of Ω, Ω/2 and related harmonics) is recovered, which is in line with the counter-intuitive trend detected through analysis of the patterning behavior. Most surprisingly a further increase in the Gershuni number (Gs = 1.00 × 10 7 , Ra ω = 1.08 × 10 6 , panel 18 (e)) causes a further simplification in the spectrum, i.e. the flow returns a single frequency response.
This peculiar evolution might be seen as reminiscent of a mechanism of the Feigenbaum type [74] occasionally mediated by intermittent events (we will come back to these important concepts in Sect. 5). Nevertheless, we could not identify a new period doubling bifurcation for the final transition from the last harmonic state to the next turbulent case (panel 18(f)); additional simulations at intermediate values of Gs and Ra ω , e.g., Gs = 3.00 × 10 7 , Ra ω =1.88 × 10 6 have simply revealed a single-frequency state akin to its predecessor for Gs= 10 7 , Ra ω =1.08 × 10 6 .

Low frequencies
While a fairly comprehensive picture has been elaborated in the previous two sections about a high frequency (Ω=O(10 4 )) and an 'intermediate' frequency (Ω=O(10 3 )), for the sake of completeness, this section continues this investigation by probing the role of a low frequency (Ω = 93.8). It can be anticipated that these cases will require a more contained discussion given their relative simplicity with respect to the more complex dynamics examined for larger values of Ω. In fact, for half of the cases simulated here (yet for Gs up to 5.00 × 10 7 ) the instantaneous field has been found to display a single possible state (in the range 10 5 <Gs<10 6 ). This class of solutions takes the form of the (aa) mode with a single cell, undergoing a periodic swap in the sense of orientation (from the clockwise orientation to the counterclockwise one, as illustrated in Fig. 19). Only a minor change in behavior can be seen in Fig. 19, panel (c), where a slight lengthening of the central roll in the x (vertical) direction is observed in the 4 th and 9 th snapshot of the sequence.
As a concluding remark for this section, we wish to highlight that recently Boaro and Lappa [75] have expanded this specific line of research (thermovibrational flow in the low-frequency regime) by addressing the companion problem in which the fluid has a viscoelastic behavior. For Pr=10, Ω = 100, comparable values of Gs, and limiting conditions for which the elasticity of the fluid is negligible (Newtonian case), they have found columnar modes of convection for Gs = 5.00 × 10 6 , 1.00 × 10 7 and 5.00 × 10 7 somehow similar to those reported here, which indicated that these are the preferred modes of convection in this range of frequencies and Pr=O(10) .

Thermal characterization: TFD distortions
After analyzing thermovibrational convection from a fine-grained mechanical level (the physicist's point of view, leading to treat the problem on the same footing as a classical study of instability or transition to chaos), this section is used to provide the reader with a more 'engineering' point of view, i.e. a coarse-grained macroscopic perspective connected to the evaluation of factors by which the thermal behavior of the system can be characterized 'globally'.
In particular, in line with the existing literature on this specific form of convection, we use the so-called thermofluid-dynamic (TFD) distortions. These can be introduced as: where T diff indicates the temperature field that would be established in the absence of any driving force (i.e. in purely quiescent conditions). Recalling Eqs. (9) and (10), Eq. (18) can also be rewritten as: where δT represents the so-called averaged distortion, i.e. and Global measures can finally be conveniently defined as .
From the literature, it is known that for thermovibrational systems, amplitudes of the fluctuating temperatures disturbances tend to decrease with frequency, and the average quantities are less dependent on the frequency so that the time-averaged distortions prevail over the fluctuating distortions at high frequencies and the fluctuating distortions prevail over the time averaged ones at low frequencies [76]. With this information in hand, it is of interest not only to explore the behavior of a fluid cell under vibrations, when the Gershuni and Rayleigh numbers are increased, but it is also of interest to study the effect of a change in the vibrational frequency.
In Fig. 21, the dashed data lines represent δT', and the solid lines represent δT. It is clear here that as the forcing frequency grows, the more the contributions of the time averaged part become predominant and vice versa, which is in agreement with the literature cited above. Interestingly, while δT' has a maximum for Ω=10 2 , δT is characterized by the presence of some peaks for Ω=10 3 (not for all values of Gs). The non monotonic behavior of the latter quantity in some circumstances can be ascribed to the complex (intermittent) dynamics, which take place for Ω=10 3 , as discussed in Sect. 4.2.
To compliment this quantitative analysis, Fig. 22 shows the timeaveraged temperature fields for the three values of Ω considered in the present work and a representative fixed value of the Gershuni number, Gs=10 6 . Although the temperature field is relatively undisturbed for the lowest frequency, as expected and in line with the trends reported in Fig, 20, for cases (b) and (c), the increase in Ω leads to an increased distortion.

Discussion and conclusions
The mechanisms through which a closed system undergoes transition to chaos can be categorized into three main types [47,70]. Firstly the Table 3 Summary of system responses for increasing values of Gs and Ra ω at Ω = 9.38 × 10 2 .

Gs
Ra ω Flow regime 1.00×10 5 1.08×10 5 Synchronous 1.70×10 6 4.47×10 5 Periodic doubling 5.83×10 6 8.27×10 5 Turbulent 6.25×10 6 8.56×10 5 Periodic doubling 1.00×10 7 1.08×10 6 Synchronous 5.00×10 7 2. Ruelle-Taken-Newhouse scenario [67,68] where a flow can be assumed to be chaotic after three incommensurate frequencies have been produced (the ensuing evolution of the system generally consisting of a corrugation of the related T 3 torus). The second known route to chaos is the so-called Feigenbaum sequence [74] where, following an increase of the control parameter, a series of period doubling bifurcations takes place. Finally, the Manneville and Pomeau [77] mechanism, where intermittent bursts occur in the time series of characteristic quantities (temperature or velocity); these bursts become more and more frequent with an increase in control parameter, resulting ultimately in a signal composed of only turbulent burst (giving rise to a fully turbulent frequency spectrum). It is known that these possible routes are not mutually exclusive, nor are they truly progressive, that is, they can mutually 'contaminate' one another.
In an attempt to fill a gap still affecting the existing literature, in the present work, the focus was on the evolution of pure thermovibrational flow. Given the known dependence of this form of thermal convection on the frequency of the applied forcing (leading to partition the space of parameters into two main regions, namely, the ranges of small and high frequencies), we have allowed such a frequency to span a relatively wide interval of orders of magnitude.
Starting from the case for which the frequency is high (a common point of origin from which many past studies have departed owing to the inherent simplifications allowed in terms of governing equations) [53], and assuming a fixed range of values of the Gershuni number (known to be the main parameter controlling the evolution of the system in the time-averaged space), the frequency has been gradually decreased with a three-fold purpose; firstly, to discern the influence of Ω on the magnitude and patterning behavior of the time-averaged velocity (especially for high values of the Gershuni number for which no results exist in the literature); secondly, to reveal the underlying textural transitions affecting the instantaneous flow field, and thirdly, to decode the typical route of evolution towards chaos.
Remarkably, it has been found that the two aforementioned regimes (generally distinguished according to the relative importance of timeaveraged and fluctuating effects) also differ in regard to the preferred sequence of events (the hierarchy of instabilities), which take place as the control parameter is increased.
For high frequencies (Ω≅ 10 4 ), the present numerical results have shown that thermovibrational convection considered in the timeaveraged space (i.e. in terms of V) is initially steady and it becomes time-periodic if the Gershuni number is increased. In this case, the typical mechanism of transition to chaos is represented by the classical Ruelle-Takens path of evolution (requiring the presence of three independent frequencies in the instantaneous velocity signal [69].) For intermediate frequencies (Ω≅ 10 3 ), however, the scenario is more complex. The ranges of existence of the time-averaged modes of convection corresponding to the quadrupolar field and the inversional symmetry pattern dramatically change. A significant modification can also be spotted in the chain of transitions leading to the development of chaos. While the hallmark of the high frequency regime is its adherence to the standard Ruelle-Takens route (which has emerged as a leading candidate to interpret most of turbulent phenomena in nature), for Ω≅ 10 3 , concurrent aspects of the Feigenbaum and Manneville and Pomeau mechanisms can be recognized. We wish to remark that, in order to eliminate the possibility that these results ensue from numerical error or solver inaccuracy, these cases were also ran using the in-house code developed by Lappa and Burel [34], and the exact same solutions obtained with OpenFOAM were found. According to both computational platforms, the uprising towards chaos as a result of the increase of the control parameter Gs (or Ra ω ) is via the Feigenbaum mechanism and the related period-doubling bifurcations; nevertheless, occasional back transitions from a chaotic condition to a sub-harmonic state can also be observed.
As a concluding remark, we wish to highlight that similar behaviors have been reported for a set of companion problems not necessarily involving the presence of an external forcing. For instance, when simulating a converging diverging channel flow, taking the Reynolds number as the control parameter, Guzman and Amon [72] observed frequency-locked periodic solutions, where the periodic behavior of the flow is restored for a particular intermediate value of Re. This was interpreted as the result of a synchronization of two oscillators spontaneously produced by the system. Similar events were also reported by Berge et al. [78].Even more relevant examples can be found in the literature about problems where a horizontal temperature gradient interacts with a constant-gravity field modulated by small harmonic oscillations. As an example, for the square cavity case, Farooq and Homsy [79,80], demonstrated that, under certain parametric conditions of finite frequency, resonant states can be induced by the interplay of convective motion of vibrational origin with the fundamental instabilities of the base flow induced by the constant gravity (similar findings are also due to Chen and Chen, and Kim et al.) [81][82][83]. In such a context, it is also worth citing Lizée and Alexander [47] who examined the case with gravity perpendicular to the applied temperature gradient and vibrations in the horizontal direction and reported on the existence of period-doubling transitions, periodic windows, strange attractors, and intermittencies as a function of the classical gravitational Rayleigh number. Keeping in mind also the arguments about the relevance of the concept of anti-resonance elaborated by Boaro and Lappa [73] in regard to the stabilization of flows which admit multiple oscillators there is no doubt that in addition to its relevance to space research and microgravity science, there are other benefits tied to a better understanding of these attractive 'vibrational' fluid systems. Looking forward, one may see turbulent thermovibrational fluid motion as a tool for finding exotic states of convection that, for now, have not been encapsulated yet in the major theories for transition to chaos.