Spatial and temporal evolution of three-dimensional thermovibrational convection in a cubic cavity with various thermal boundary conditions

Thermovibrational ﬂow in a differentially heated cubic cavity with vibrations applied in a direction parallel to the imposed temperature gradient is investigated by solving numerically the governing equations for mass, momentum, and energy in their original nonlinear form. A parametric analysis is conducted through the stepwise examination of the following degrees of freedom: magnitude of the Rayleigh number and the thermal behavior of the sidewalls. A complete characterization of the emerging time-varying convective structures is attempted in terms of spatial symmetries broken or retained, related temporal evolution, and global parameters, such as the Nusselt number. It is shown that the intrinsically three-dimensional nature of the problem and its sensitivity to the thermal boundary conditions can have a remarkable inﬂuence on the multiplicity of emerging solutions and the system temporal response


I. INTRODUCTION
The investigation into natural convection has been under way for centuries.Since the works of Hadley, 1 providing the first theory explaining the behavior of convection currents on our planet or trade winds, and the seminal experiments with shallow liquid layers placed on a heated wall conducted by Rayleigh, 2,3 the study of this specific kind of fluid motion has become wide-spread due to its omnipresence in nature and relevance to industrial processes.As a result of the aforementioned landmark investigations, the general field of buoyancy convection is nowadays overarched by two archetype configurations: the Hadley flow, where a fluid cavity is heated from the side and subjected to vertical steady gravity (often simply referred to as the laterally heated-cavity problem) and Rayleigh-B enard (RB) convection where fluid is heated from below (yet in the presence of a steady vertical gravity field).
In addition to these two forms of classic thermo-gravitational convection, another variant (discovered much more recently) can be obtained if in place of changing the relative direction of the temperature gradient and gravity, the latter is replaced with a time-varying acceleration.This specific kind of fluid motion is generally known as thermovibrational convection.Notably, when gravity is replaced with (or superimposed onto) vibrations, in addition to the standard practice of varying the strength of the resulting flow by increasing or decreasing the temperature gradient across the system, two additional degrees of freedom are unlocked.As the vibrations are defined as a sinusoidal displacement in time, the flow they can produce can be modulated by varying their amplitude and/or their frequency.
Initial studies on thermovibrational convection date back to the 1970s (Gershuni and co-workers [4][5][6][7][8][9][10] ) when the concept of weightlessness started to attract interest due to the availability of relevant spacecraft for the execution of experiments in microgravity (sounding rockets and later the space shuttle).Nowadays, studying the effects of vibrations on fluid systems is of paramount importance as they are representative of the influence that various inertial disturbances (g-jitters) may have on available space platforms, such as the International Space Station (ISS).These residual gravitational disturbances are of different types.In particular, according to the classification elaborated by Nelson 11 (who considered oscillatory disturbances that are essentially harmonic and periodic in nature), these can be quasi-steady (i.e., aerodynamic drag solar radiation pressure, etc.), transient (vehicle maneuvers, crew motions, etc.), and oscillatory (structural vibrations, operation of machinery, etc.).The ISS is equipped with various systems to mitigate these perturbations as they can be disruptive to fluid-dynamic experiments (for a relevant analysis of the tolerability of fluid experiments to them, the interested reader may consider, e.g., Refs.12-17).
Most remarkably, following these attempts, the field of thermovibrational convection has evolved from its initial heartland of spaceexperiment-supporting research into a completely independent discipline attempting to discern how the nonlinearities present in fluid systems conspire to form organized spatial patterns.While the g-jitters that occur on space platforms are typically high frequency and low amplitude vibrations and therefore only occupy a small region of the possible space of parameters, recent studies on these subjects have revealed that when the entire spectrum of vibrations is considered, many new fluid responses and related instabilities become possible, which (due to their peculiar properties) have proven to be irresistible to theoretical physicists.
Works treating this particular field of study can yet be separated into two main classes akin to the two categories previously mentioned (Hadley flow and RB convection) where vibrations are either perpendicular or parallel to the temperature gradient across the cavity.6][27] As for what concerns the parallel configuration, Hirata et al., 28 and more recently, Crewdson and Lappa 29 have provided a comprehensive map of possible system states and patterning behaviors for a wide range of vibrational modes and Rayleigh numbers.
Unfortunately, however, these studies have been limited to twodimensional (2D) geometries, thereby filtering out the inherent complexities typical of effective three-dimensional (3D) flows.In order to fill this gap, in the present work the set of predictive links between flow properties and related influential factors is expanded through consideration of a real (3D) cavity.Given the lack of equivalent results in the literature, in particular, a cubic enclosure is considered.This specific choice is motivated by the availability of a significant amount of existing data for the canonical case of standard Rayleigh-B enard convection (i.e., cubic enclosures uniformly heated from below and cooled from above with vertical steady gravity), which can be considered for comparison and/or as a guide to interpret the still completely unknown behavior of thermovibrational flow in similar conditions.
Along these lines, we wish to recall the earlier investigation by Pallare `s et al., who addressed the RB problem both numerically 30 and experimentally. 31In their studies, seven fundamental modes of convection were identified for fluid motion driven by gravity in parallelepipedal cavities heated from below (these including single-roll, two-roll and other toroidal roll-type states; we will provide a more detailed description of these structures at a later stage in the present work).The stability of these solutions was found to depend on the Rayleigh number (Ra), Prandtl number (Pr) of the fluid, and the aspect ratio of the cavity.Puigjaner et al. 32,33 conducted a stability analysis of the same problem and found that by increasing the strength of convection (i.e., Ra), a series of possible steady, stable (ST), and unstable patterns emerge resulting in a complex bifurcation diagram.Despite some differences in the related findings, the greatest merit of these valuable efforts resides in the provided evidence that an apparently innocuous problem, such as thermogravitational convection in a 3D cubic cavity, can offer a relevant playground for the emergence of complex behaviors and the analysis of the underpinning nonlinearities.
Lappa 34 argued that the discrepancies affecting Pallare `s et al. 31 and Puigjaner et al. 32,33 should be justified when considering that the transition from one solution to another, upon an increase in the Rayleigh number, is dependent on the presence of thermal and/or momentum boundary layers and the relative ratio of their thickness, which in turn changes with the value of Pr and the specific thermal behavior of the sidewalls.The effect of such conditions on the emerging flow structures has been investigated experimentally by Pallare `s et al., 31 for Pr ¼ 130 over a range of moderate Rayleigh numbers (Ra < 8 Â 10 4 ).Both adiabatic and conducting wall conditions were found to yield similar results in terms of flow patterns and transitions up to Ra % 5 Â 10 3 .Beyond this threshold, however, the flow structures were observed to differ depending on the considered thermal conditions and to occur at different values of the Rayleigh number.In this regard, we wish also to cite the later study by Puigjaner et al., 35 who considered fully conducting walls for a wide range of Ra and Pr, and provided relevant arguments about the influence of the thermal conditions imposed on the sidewalls.
Motivated by this observational tide, we therefore concentrate on thermovibrational convection in a cubic cavity driven by shaking imposed in a direction parallel to the temperature gradient and assume different types of thermal boundary conditions for the sidewalls.Specifically, in addition to the canonical configurations where these are either completely adiabatic or conducting, a hybrid situation is also investigated where, while two opposing walls are adiabatic, the other couple maintains a linear temperature profile.The ensuing numerical results are presented in terms of flow structures, patterning behavior and the (often not obvious) relationship that is established between the forcing (the vibrations) and the temporal response of the flow.

II. MATHEMATICAL MODEL
The simple 3D cubic cavity considered in this study assuming microgravity conditions is shown in Fig. 1.The bottom wall is cold and the top wall is hot.As anticipated in Sec.I, the vibrations are applied parallel to the temperature gradient and three thermal wall conditions are investigated [Figs.1(a)-1(c)].

A. Governing equations and their non-dimensional form
The vibrations are modeled mathematically as a sinusoidal displacement in time, characterized by the amplitude of the displacement b (m) and the angular frequency of the vibrations x (rad/s), where x ¼ 2pf and f is the frequency in Hz.Mathematically, this gives where n is the direction vector of the vibrations (unit vector).The time-varying acceleration (to be accounted for in the production term of the momentum balance equation) is then obtained by taking the second derivative of Eq. (1), g ðtÞ ¼ g x sin ðxtÞn ; where Integrating g(t) over one vibrational period (2p/x), one gets This indicates that the time-averaged value of the acceleration over one vibration period is equal to 0. The balance equations for mass, momentum, and energy under the assumption of incompressible flow read where V, p, and T are the problem primitive variables, i.e., velocity, pressure, and temperature, respectively.Moreover, q is the fluid density, l is the fluid dynamic viscosity, and a is the fluid thermal diffusivity.Using the Boussinesq approximation, the density is considered everywhere constant with the exception of the momentum production (buoyancy) term at the right-hand side of Eq. ( 5), which is expanded as where T o is a reference temperature, q 0 is the density for T ¼ T o , and b T is the so-called thermal expansion coefficient by which variations in the fluid density can be directly linked (via a linear relationship) to the gradients of temperature in the fluid.By substituting this term into the momentum equation, assuming that the individual term q 0 g is absorbed into the modified pressure term and scaling all lengths by the size of the cavity (L), the velocity by (a/L) the time by (L 2 /a) and the pressure by (q 0 a 2 /L 2 ), all these equations can be cast in compact (non-dimensional) form as where Pr is the well-known Prandtl number, ( being the fluid kinematic viscosity given by the ratio of dynamic viscosity and fluid density, i.e., ¼l/q 0 ).Moreover, Ra x appearing in the buoyancy term is the vibrational Rayleigh number, analogue to the classical Rayleigh number used in standard gravitational convection problems (Ra), namely, Finally, X is the non-dimensional angular frequency of the vibrations defined as where P is the non-dimensional period of vibrations.

B. Boundary and initial conditions
At t ¼ 0, the velocity field across the cavity for all cases is 0, Moreover, a linear temperature profile is assumed along y (this reflects the widespread practice of performing experiments in space on fluids subjected to vibrations after having "thermalized" the fluid container, i.e., after having established a diffusive distribution of temperature in order to filter out transient thermal effects).The boundary conditions for t > 0 are reported in the following.
No-slip conditions are applied for all the solid walls everywhere.Moreover, different thermal boundary conditions are considered depending on the chosen setup (as shown in Fig. 1).In particular, for the case illustrated in Fig. 1(a) the sidewalls are adiabatic @T @x ¼ 0 for z ¼ 0; z ¼ 1 and t > 0; (15a) @T @z ¼ 0 for x ¼ 0; x ¼ 1 and t > 0: (15b) For the case illustrated in Fig. 1(b), the sidewalls in the z-plane are perfectly conducting and the walls in the x-plane are adiabatic, @T @x ¼ 0 for z ¼ 0; z ¼ 1 and t > 0; (16a) For the case illustrated in Fig. 1(c), all the sidewalls are perfectly conducting, i.e., T ¼ y for x ¼ 0; x ¼ 1 and z ¼ 0; z ¼ 1 and t > 0: (17)   Finally, different constant temperatures are set for all three cases at y ¼ 0 and y ¼ 1, T ¼ 0 for y ¼ 0; 0 x 1; 0 z 1 and

III. NUMERICAL METHOD
The numerical simulations presented in this work have been produced by solving the balance equations presented in Sec.II [Eqs.( 8)- (10)], using a finite volume method for incompressible flow.In particular, the computational platform OpenFOAM has been used.This relies on the well-known PISO algorithm (pressure implicit split operator) to ensure adequate pressure-velocity coupling.The details of this time marching procedure can be found in a number of works previously published by the present authors (see, e.g., Ref. 29).However, for completeness, the basic principles of the method are also briefly reported here.

A. Pressure-based approach
The PISO algorithm pertains to a class of projection methods, where in the first instance the pressure (p) is artificially removed from the momentum equation giving rise to the following equation: where only the velocity field requires solving for.This velocity field is, however, only provisional as it does not account for pressure and is simply used to initiate the time marching procedure.It also does not satisfy the incompressibility constraint set by the continuity equation [Eq.( 8)].To overcome this, the Hodge decomposition theorem is used to split V Ã into two contributions.This theorem states that any vector field can be decomposed into a divergence-free contribution and the gradient of a scalar potential, giving where C is a constant.As p is still unknown and assuming C ¼ Dt, Eq. ( 20) can then be substituted into the continuity equation ( 8) giving This equation can be used after the integration of Eq. ( 19) to calculate the pressure.If the effective physical velocity is imposed on the boundary of the domain hosting the fluid, in particular, homogeneous boundary conditions can be considered for it, i.e., @p/@n ¼ 0 (where in this case n denotes the direction perpendicular to the boundary).Once the pressure is known, a "physical" divergence free velocity is finally calculated using Eq.(20).As for what concerns the numerical schemes employed, we have exploited upwind and central differencing schemes in space for the convective and diffusion terms, respectively.Moreover, a multicore processing approach has been implemented where the domain is decomposed into four sub-domains (by cutting it with transverse perpendicular planes located at x ¼ 0.5 and z ¼ 0.5, respectively).A fixed time step has been used for the time integration, satisfying the condition Dt ¼ 2p=N p x where N p is the number of time steps over the vibrational period (P ¼ 2p=x).By setting N p ¼ 3257, the twofold objective of ensuring reasonable solution accuracy and algorithm numerical stability has been met.

B. Validation and grid refinement
A cross-validation approach has been adopted in the present work because limited literature or benchmark solutions exist for thermovibrational convection (especially for the case where vibrations are parallel to the temperature gradient; an exception being the authors previous study, 29 which however was validated for a different value of Pr).Given this drawback, most conveniently, the ability of OpenFOAM to capture properly the dynamics of interest has been verified through comparison with the results obtained using an in house explicit finite volume solver. 34For simplicity and to save computational time, this assessment has been carried out for a 2D square cavity with conducting sidewalls.The average Nusselt number across the hot wall and both the vertical and horizontal velocity components at the point (0.75,0.75) of the cavity have been compared.
As quantitatively substantiated by Table I and the velocity signals reported in Figs. 2 and 3, a satisfactory agreement has been found (we wish to highlight that the small discrepancy affecting the velocity signals determined with the two different computational platforms must be ascribed to the different location of the primitive variables on the computational grid, this being based on a collocated and a staggered grid approach in OpenFOAM and the explicit finite volume solver by Ref. 34, respectively).
Thereafter, a mesh refinement study has been carried out on the 3D cavity with adiabatic sidewalls for Ra x ¼ O (10 6 ).In order to make the overall process more efficient, some theoretical criteria have been applied in the attempt to estimate "a priori" the required number of grid points.First, given the tendency of thermal convection to display turbulent behavior when the Rayleigh number attains high values, the so-called Kolmogorov length scale has been evaluated as a possible In such a context, it is worth recalling that this length scale provides a restrictive condition in terms of non-dimensional size of the mesh (Dx c ), i.e., it is expected to satisfy the condition Dx c < f Ra .This constraint can be easily used to determine accordingly the number of divisions required across the domain as.N div1 ¼ 1=f Ra .Assuming the worst condition, i.e., the maximum value of the vibrational Rayleigh number being considered in the present work, Ra x ¼ 8.34 Â 10 5 , N div1 ¼ 1=ð1:61 Â 10 À2 Þ ' 62. From this, we have deemed it necessary to perform a mesh refinement study allowing the number of divisions to range in a certain neighborhood of this value.Moreover, since it is known that for these high values of Ra and a relatively high value of Pr, thermal boundary layers are prone to develop across the top and bottom wall of the cavity, the number of cells required in this boundary layer has also been taken into account following the criteria provided by Shishkina et al., 38 N BL ffi 0:35Ra 0:15 : A correlation about the thickness of the thermal boundary layer can be found in the earlier study by Russo and Napolitano, 39 From these two relationships, we have calculated: N BL ¼ 2.70 and d thBL ¼ 3.31 Â 10 À2 .Using these two values, for a uniform grid, the number of divisions has finally be determined as which gives N div2 ¼ 81.75.With these theoretical requirements in hand, a mesh refinement study has been carried out.Including first a mesh containing 60 3 elements and followed by 80 3 , 100 3 , and so on.The ensuing results shown in Fig. 4 have been taken when the flow is the most disturbed.It can be noted that the data for the 60 3 and 80 3 are similar; however, as in these cases the theoretical requirement outlined by N div2 ¼ 81.75 is only just satisfied, a further jump of 20 divisions has been considered.As from 80 3 to 100 3 a significant variation in value can still be observed, a further jump to 120 3 has been implemented.The results obtained with the 100 3 and 120 3 meshes are indecipherable, and therefore, the 100 3 mesh has finally been chosen as the required mesh for Ra x ¼ 8.34 Â 10 5 .

IV. RESULTS
In this section, the outcomes of six individual simulations are presented.These cases encompass the three types of heating configurations defined in Sec.II B for two different orders of magnitude of Ra x (namely, Ra x ¼ 8.34 Â 10 4 and Ra x ¼ 8.34 Â 10 5 ).Before starting to deal with such findings, we wish to recall that, as shown by Hirata et al., 28  constraint of two-dimensionality, a rich map of solutions is possible in systems where the vibrations are set parallel to the temperature gradient when both the vibrational Rayleigh number and frequency of vibration are varied.By allowing these systems to develop in a realistic 3D space, an even larger set of states must therefore be expected.In turn, these can differ in regard to the instantaneous spatial organization (the flow structure and associated symmetries) and the related evolution in time (the system "temporal" response).Given the inherent complexity of the problem, in the following a peculiar approach is implemented where an attempt is made to treat these (spatial and temporal) aspects in a separated manner (however, still creating the relevant links as necessary).
The cases are setup as per Fig. 1, and the fluid parameters are shown in Table II.
Given the otherwise intractable scale of the problem, without loss of generality, the angular frequency of the vibrations is fixed to X ¼ 50, a value for which (as illustrated in Sec.IV A) the flow is expected to develop a remarkable degree of unsteadiness.

A. Fluid response and velocity signals
Before entering into a discussion regarding the threedimensional textural transitions affecting the patterning behavior, in Sec.IV A the velocity signals obtained from probes (similar to those exploited for the validation study) are used to identify the regime embodied by the flow.Along these lines, it should be pointed out that four possible solutions or regimes have previously been recognized in the 2D study by Hirata et al.: 28 synchronous (SY), subharmonic (SU), non-periodic (NP), and stable (ST).Two additional solutions were reported in the later analysis by Crewdson and Lappa, 29 namely, the synchronous and periodic or synchronous and non-periodic (SY-P or SY-NP) solutions.
Two of these possible behaviors can yet be recognized in the results obtained for the 3D configuration examined in the present work.These include the SY-P and SU regimes.For the former (SY-P), the flow repeats itself periodically.For the latter (the SU case), the frequency of repetition of the fluid behavior is halved with respect to the forcing frequency, resulting in a period twice as long as the forcing period.
In particular, Fig. 5 shows the velocity signals for the lower value of Ra x considered here (Ra x ¼ 8.34 Â 10 4 ), where the SY-P regime is evident for all the three variants of thermal boundary conditions defined in Sec.II B.
Additional insights can be gathered from Fig. 6, where it can be seen that the signals appear to be of the type SY-NP (displaying turbulent bursts) as the Rayleigh number is increased by one order of magnitude.A closer inspection of these signals further reveals that a SU regime is enabled for case (b) (indeed, the periodicity in shape of the signal is repeated every two periods in this case).
On the basis of this initial assessment (relying solely on arguments based on the time response of the flow), it may therefore be concluded that when 3D configurations are considered, a change in the thermal boundary conditions can cause a remarkable variation in the temporal behavior of the flow.As illustrated in detail in Secs.IV B and IV C, such responses are intimately coupled with a variety of textural transitions in the flow structure, which deserve their own treatment.

B. Flow structure characterization for small Rayleigh numbers
In order to support the reader's understanding of such dynamics (before diving into a purely spatial characterization of them), it is worth recalling some key aspects of thermal convection which will prove very useful later for their interpretation.In particular, it is instructive to recall that this kind of fluid motion is governed by the same equations everywhere in space, yet it takes a form that has periodic spatial variations, with "nodes" (velocity or temperature extremes) positioned at given points.Physicists say that it "spontaneously breaks space-translation symmetry."As already explained to a certain extent in Sec.I, this useful concept can be used as a tool to characterize or categorize the flows in certain universality classes, which transcend the specific nature of the considered flow (be it driven by steady gravity or vibrations).Section IV B 1 is, therefore,  dedicated to the provision of some necessary propaedeutical arguments along these lines.

Textural transitions and patterning behaviors
Having completed a sketch of the observed temporal response of the flow in Sec.IV A, the focus is now shifted to the specific spatial behavior of the fluid over one vibrational period.Along these lines, starting with the case where Ra x ¼ 8.34 Â 10 4 and all four sidewalls are set as adiabatic, Fig. 7 shows the related evolution of the fluid streamlines.As the reader will immediately realize by inspecting this figure, most conveniently, we have reported the sequence of snapshots of the flow structure in conjunction with the signal already shown in Fig. 5.This is instrumental in making evident that the degree of complexity displayed by the system strongly depends on the considered specific sub-region of the period.In particular, a relatively wide subinterval exists where convection can be considered relatively weak (almost quiescent situation, hereafter referred to as "quasi-stationary state"), whereas its amplitude greatly grows as conditions are examined that correspond (or are located in proximity to) to the signal peak.
Accordingly, we split the analysis into two parts, the first part, being the time over which the fluid is quasi-stationary (where the flow adopts a resting configuration and the velocity of the fluid is minimal) and second, the time region where the aforementioned convective pulse occurs.Looking at Fig. 7, the first and last frames show that the resting configuration of the flow is essentially a singular toroidal roll.
During the convective phase, the toroidal roll symmetry is maintained; however, the roll becomes more compact at first and then undergoes a series of minor textural transitions before regaining its resting configuration.In particular, a more detailed impression of the streamline behavior at t ¼ 0.3P as well as the distribution of the related velocity components U x , U z and U y in the xz plane (at y ¼ 0.5) can be gathered from Fig. 8.
Toward the end to identify universality classes in these behaviors, reference could be made to meaningful classifications introduced in the past for classical RB convection.Relevant examples of this modus operandi can be found, e.g., in the work by Mizushima 40 and Mizushima and Adachi 41 for the case of 2D RB flow.
When moving to the 3D case (i.e., the cubic cavity shown in Fig. 1), however, the simple criteria valid for 2D configurations become rather inadequate (only partially able to account for certain properties of the flow).A more exhaustive characterization approach may, therefore, be based on the (seven) fundamental modes that Pallare `s et al. 30 originally identified for standard RB convection in a cubic enclosure in the range 3.5 Â 10 3 < Ra < 6 Â 10 4 for Pr ¼ 0.71, 10, and 130.These solutions are denoted in the following with Sn with n being an integer ranging between 1 and 7; see Table III and Fig. 9.
Additional structures were identified by Puigjaner et al.; 32 however, a description of their geometric qualities is not provided here (these included the modes S9 and S11-S15).
As the reader will realize by inspecting Table III and Fig. 9, these solutions essentially reflect the symmetries of the group D 4h ¼ Z 2 Â D 4 .In particular, the sub-group D4 includes the symmetries of  a regular polygon with four vertices.In turn, these consist of the mirror reflections with respect to the middle (x ¼ 1/2 and z ¼ 1/2) and diagonal (x ¼ z and z ¼ 1Àx) vertical planes and related combinations.
In the framework of this sub-group, the solution corresponding to the toroidal roll (dominant in the quasi-stationary regime) might be, therefore, considered as a fundamental mode S4 (Pallare `s et al., 30 see Table III).The symmetry with respect to the horizontal mid-plane (y ¼ 0.5), however, must be also considered.Unlike the other symmetries pertaining to the above-mentioned dihedral sub-group D 4 , which are also applicable to the temperature field, this symmetry obviously applies only to the velocity field (it is equivalent to rotations of the velocity field of an angle p around one of the x or z horizontal directions in the y ¼ 1/2 plane; the reader may consider Puigjaner et al., 33 for an exhaustive mathematical description of all these groups of symmetry, which is not reported here for the sake of brevity).
Looking forward to the next cases, we omit the representation of the velocity signal as, from Fig. 5, we know that the three cases display a high degree of similarity from a temporal point of view.Figure 10 illustrates the corresponding evolution of convective modes for the other two boundary condition configurations (still for Ra x ¼ 8.34 Â 10 4 ).
In particular, in Fig. 10(a), a planar configuration along the z axis of four separate rolls can be recognized.This quadrupolar roll configuration is a well-known solution in 2D thermovibrational studies, reported early on by Gershuni et al. 6 A somehow similar steady state has been found by Pallare `s et al., 30,31 (S5), described as a four roll structure, where however each roll axis is perpendicular to one sidewall.In the present case [Fig.10(a)], the axes of rotation of the rolls are perpendicular only to the two conducting sidewalls, thus resulting in a planar structure rather than the S5 mode where symmetry was found about the diagonal x ¼ z.A graphical representation of the resting configuration shown in Fig. 10(a), akin to the classification style of Pallare `s et al. 31 is provided in Fig. 11.
The symmetry with respect to the y ¼ 1/2 midplane is embodied here where an even number of rolls is found along the x and y axis, with no rolls found along the z axis due to the planar (almost 2D) nature of the flow.
The advent of the convective phase, for the case shown in Fig. 10(a), brings about a transition from the four-roll to a two-roll planar convective configuration followed by the planar structure disappearing past t ¼ 0.3P.This next convective mode (shown in more detail in Fig. 12) is the result of the fluid motion attempting to overcome the constraining boundary conditions and form a toroidal structure, similar to that seen for the fully adiabatic situation (however, bound by the sidewalls a fully toroidal structure is not achieved).Finally, the resting configuration is reestablished.
The next figure of the sequence (Fig. 13) refers to the case where all four sidewalls are conducting.As a fleeting glimpse into this figure [and Fig. 10(b)] would confirm, during the resting phase, the flow adopts a hybrid configuration of the fully adiabatic case and the half adiabatic, half conducting case.The toroidal S4 structure seen in Fig. 9 is still visible, however, the cavity hosts two torii instead of one.In this case, the symmetry with respect to the y ¼ 1/2 plane can be therefore   recognized if the flow is observed from the front and side view of the cavity.For what concerns the transitional behavior during the convective phase, the stages of evolution are similar to those found for the fully adiabatic case; however, past t ¼ 0.3P, the dominant convective mode is compressed at the lower half of the cavity preceding the advent of small rolls eventually responsible for the two separate torii visible in the resting configuration.
On the basis of this initial set of results, the following conclusions can, therefore, be drawn.When comparing the effects of the three types of boundary conditions on the flow, it becomes evident that the fully adiabatic and fully conducting cases bare much similarity as symmetry about the x and z axis at the zx mid-plane is maintained through both the resting and the convective phases.However, an additional transitional structure is observed for the fully conducting sidewall case.By contrast, the case with two adiabatic and two conducting sidewalls is different as the flow structure is essentially planar in the quasi-stationary regime.
Another striking analogy applies to the configuration with the conducting walls.For both the hybrid case and the fully conducting sidewall case, an antisymmetry emerges along the y axis.This leads to the classic quadrupolar field in the case with hybrid thermal boundary conditions (displaying a planar four roll structure) and produces a dual-toroidal structure in the fully conducting situation.
However, as the reader will easily realize by taking a closer look at Figs. 7, 8, 10, 12, and 13, due to the complexity of the structures observed during the convective phase (when the fluid is not at rest), the application of earlier classifications, such as those developed by Pallare `s et al., 30,31 is not always possible.
For this reason, in Secs.IV B 2 and IV C, a further level of analysis is implemented through the consideration of specific aspects emerging from observation of the profiles of the different velocity components and/or the related distributions (maps) in certain planes as a function of time.
In particular, Fig. 14 shows the lines over which the velocity profiles are taken.The x, y, and z velocity components (U x , U y , and U z ) are taken over both (dashed) centerlines for all three cases.

Velocity profiles and emerging symmetries
Starting with the fully adiabatic case, Fig. 15 depicts the velocity profiles taken at intervals of 0.1P over one vibrational period.The highest velocity magnitude occurs at t ¼ 0.3P for both the vertical and horizontal components.Figure 15(a) shows a positive U x over the first half of the x axis and a negative U x over the second half of the x axis, indicating that the center of the cavity acts as an attractor for the duration of the entire period.The profile displayed in Fig. 15(b) exhibits a positive vertical velocity U y at the center of the cavity indicating an upward motion of the fluid.This is in agreement with the arguments provided earlier about the flow observed to be circulating up and outward from the center of the cavity.
The U z velocity profiles are not provided, and neither is the U x velocity along the z axis, owing to the fact the flow adopts the toroidal structure, and therefore such information would be redundant (U x along the x axis ¼ U z along the z axis, and U x ¼ 0 along the z axis ¼ 0 and U z along the x axis ¼ 0 and finally U y along the x axis ¼ U y along the z axis; Fig. 16 shows the velocity fields confirming that these equalities hold for y ¼ 0.5).
For the situation where two sidewalls are conducting and the front and back walls are adiabatic, we look first at the horizontal velocity components (U x ) and (U z ).Unlike the fully adiabatic case, this case does not display perfect agreement between the U x and U z velocity profiles.This is seen in Fig. 17 and confirmed in Figs.19(a) and 19(b).The vertical velocity profile, however, can be seen to vary significantly depending on the sampling axis, as shown in Fig. 18.Having found many similarities between the fully adiabatic and fully conducting case during the convective pulse stage, related results are not described in this section for brevity.

C. High Rayleigh number
Looking now at the second value of the vibrational Rayleigh number considered in the present study (Ra x ¼ 8.34 Â 10 5 ), it is worth starting from the simple remark that similarly to the previous case, the flow embodies either a quasi-stationary state or a convective burst (this being illustrated in Fig. 20).The convective burst, however, displays a much more complex behavior as the fluid is more disturbed during this part of the period in comparison with the equivalent behavior seen for the lower value of Ra x .If the Rayleigh number is increased by one order of magnitude, a myriad of solutions emerge and disappear.
Here, we look first at the resting configurations for all three boundary conditions.In particular, as qualitatively substantiated by Fig. 21, the resting configurations for the case Ra x ¼ 8.34 Â 10 5 with adiabatic sidewalls are identical to the equivalent configurations identified in Sec.IV B for Ra x ¼ 8.34 Â 10 4 .
Figure 22 naturally complements Fig. 20 by revealing the evolution of the system over one vibrational period and the related multiplicity of solutions excited during the convective stage (burst).Looking first at the panel t ¼ 0, the toroidal structure is accompanied by a quasistationary fluid and a linear temperature distribution along the y axis.
As evident in the close-up in Fig. 23, although the fluid pattern displays symmetry about the x and z axes at the zx mid-planes, at t ¼ 0.1P the convective burst is enabled and the application of the concepts developed by Pallare `s et al. 30,31 becomes rather challenging.
Moving through the period, the symmetry of the system and the number of rolls and their orientation seem to vary in a relatively random way.The temporal behavior of the flow is reported in Fig. 24, where the velocity profiles of all three cases are shown.The velocities are taken over a line cutting through the center of the cavity as per Fig. 14.
Looking first at the vertical velocity component (U y ) taken parallel to the x axis, symmetry is visible with respect to the center of the cavity (x ¼ 0.5), where a peak in velocity is seen for all three cases, excepting at t ¼ 0.3P for the fully adiabatic wall case.This is also observed when looking at vertical velocity component taken parallel to the z axis, reported in Fig. 25(a).
With exception of this momentary break in symmetry, the vertical velocity profiles are inherently symmetrical and remain so at least until t ¼ 0.4P.After this point, the velocity of the fluid is low compared to the first instances of the convective burst (t ¼ 0.P to 0.4P).Even though the velocity is low, many textural transitions are observed (as evident in Fig. 22) while the fluid settles from the convective state to the quasi-stationary state.As for the horizontal velocity components, the temporary asymmetry observed in the vertical velocity component in the fully adiabatic case is widespread through the whole period indicating that for the case where all walls are adiabatic, the fluid is more prone to breaks in symmetry than for the situation when conducting sidewalls are considered.Looking at Figs. 26(b), 26(c), 27(b), and 27(c), symmetry is apparent about the center of the sampling line.
For the sake of completeness, following other authors, 42,43 the spatiotemporal maps providing a visualization of the symmetry (or asymmetry) embodied by the flow over time are also included (Figs. 28-30).
The asymmetry embodied by the fully adiabatic case can be recognized in Figs.28(a), 28(d), 29(a), 29(d), and 30(d).However, a near symmetrical pattern is maintained for the vertical velocity component along the x axis for the fully adiabatic case [Fig.30(a)].
Moving on to the next two cases, the half conducting half adiabatic and the fully conducting configuration, it is worth noting that both situations present a high level of symmetry for the horizontal velocity components, especially at the point where the fluid is most disturbed [as the reader will realize by inspecting Figs. , is maintained over the whole period.Finally, an additional relevant remark can be made regarding the periodicity of the hybrid case.For all maps pertaining to this case, the period over which the fluid repeats itself is twice that of the forcing period.This is in agreement with the velocity signals provided in Sec.IV A.

D. Thermal response
This section is finally devoted to an analysis of the thermal response of the system.Along these lines, Fig. 31 shows the ratio of heat transfer due to convection over that of conduction along a given boundary, i.e., the Nusselt number (Nu) for the hot wall, where Nuðx; tÞ ¼ ð 1 0 @Tðx; z; tÞ @y dz and Nu ¼ Similarly to the velocity signals, the Nu signal associated with the lower value of Ra x is of the type SY-P, whereas the responses associated with the higher value of Ra x, although still technically classifiable as SY-P, border on the verge of the SY-NP regime where a turbulent burst appears after the convective peak.It must be noted that the behavior related to the heat transfer of the system is quasi-identical when comparing the different sidewall conditions, especially when the value of Ra x is increased.
Temperature maps are also included here for all considered cases.Focusing first on Fig. 32 (Ra x ¼ 8.34 Â 10 4 ), perfect agreement is observed between the thermal behavior of the flow along the x and z axes for the situation where all four walls are either adiabatic or conducting; this is, however, not the case when a higher value of Ra x is examined (Fig. 33).
Indeed, a slight deviance is observed in the temperature profile across the z axis in this case, also detectable in the velocity maps provided in Figs.28-30.

V. CONCLUSION
Although numerous investigations have been carried out in the general area of thermal convection, an insightful and complete    understanding of the properties of the specific variant driven by timevarying accelerations has hitherto been unclear.This study is a contribution to improve the present unsatisfactory situation, especially for what concerns the poorly considered situation in which vibrations are parallel to the imposed temperature gradient.The strategy undertaken in earlier work of the present authors based on the numerical solution of the two-dimensional Navier-Stokes and energy equations has been further pursued by allowing the flow to develop in a realistic 3D physical domain for which the problem of pattern selection has long been a theoretical puzzle even for the canonical case of standard steady RB convection.

Physics of Fluids
The simulations have shown that an increase in the system (spatial) dimensionality has a dramatic influence on the richness of the fundamental modes of convection that can be excited.These can be partially grouped in different categories according to some existing classifications based on various symmetries that are broken or retained and the number of convective structures present at the same time in the physical domain.Along these lines, we have referred directly to the existing literature for RB convection in cubic cavities just to encourage readers to follow up various details and recognize the analogies and differences with this parent category of thermal flows.Given the intrinsically time-varying nature of thermovibrational convection in systems where the vibrations and the temperature gradient are concurrent, many of these convective states can be produced for a fixed value of the vibrational Rayleigh number and given thermal boundary conditions.Although two well-defined convective stages can always be identified in the period of vibrations (one corresponding to an almost quiescent quasi-stationary state, and another where a convective pulse occurs), however, the enabled modes are not mutually exclusive, nor are they truly progressive.Moreover, their multiplicity tends to be enhanced as the vibrational Rayleigh number is increased and the convective pulse is turned into a turbulent burst.
The numerical simulations have also revealed that, despite this multiplicity, some control on the morphology of the emerging convective structures can be exerted by forcing the system a priori to break or adhere to some spatial symmetries by imposing non-uniform thermal boundary conditions along the sidewalls.An ordered combination of adiabatic and conducting walls can indeed limit the ability of the flow to produce toroidal states in favor of more two-dimensional solutions.The intentional use of hybrid thermal boundary conditions can also be instrumental in inducing changes in the temporal response of these systems, causing a shift from synchronous (flow oscillating at the same frequency of the forcing) to sub-harmonic (period doubling bifurcation) behaviors or vice versa.Moreover, we have shown that cases containing a pair or two pairs of conducting sidewalls are more stable than the configuration with adiabatic sidewalls.
An exciting prospect for the future is to enrich this problem with the addition of solid particles, thereby giving rise to a new line of inquiry running in parallel with that where the interplay of thermovibrational effects and particle inertial effects in cubic cavities with vibrations perpendicular to the temperature gradient has been found to support fascinating particle self-organization phenomena (formation of highly ordered, high-resolution structures with the morphology of quadric surfaces [44][45][46][47][48] ).

FIG. 1 .
FIG. 1. Sketch of the considered geometry and related thermal boundary conditions: (a) all lateral walls are perfectly adiabatic, (b) two lateral walls are perfectly conducting, and the front and back walls are adiabatic, and (c) all sidewalls are perfectly conducting.

FIG. 4 .
FIG. 2. Vertical velocity signals (U y ) for the case Ra x ¼ 4 Â 10 4 , X ¼ 5 and Pr ¼ 7 for OpenFOAM and the in-house code.A grid size of 100 Â 100 is used.

FIG. 5 . 6 V
FIG. 5. Vertical velocity components (U y ) recorded at the center of the cavity (0.5, 0.5, 0.5) for (a) the fully adiabatic case, (b) the case where the front and back of the cavity are adiabatic and the two sidewalls are conducting, and (c) the case where all sidewalls are conducting for Ra x ¼ 8.34 Â 10 4 .

FIG. 6 .FIG. 7 .
FIG. 6. Vertical velocity components recorded at the center of the cavity (0.5, 0.5, 0.5) for (a) the fully adiabatic case, (b) the case where two sidewalls are adiabatic and two are conducting, and (c) the case where all sidewalls are conducting and Ra x ¼ 8.34 Â 10 5 .

FIG. 8 .
FIG. 8. Snapshots of flow streamlines (top) and velocity components across the mid-plane of the cavity (bottom) for the case where all sidewalls are adiabatic and Ra x ¼ 8.34 Â 10 4 at t ¼ 0.3P.TABLE III.Description of flow structures identified by Pallare `s et al., 30,31 and Puigjaner et al. 32 Structure names Pallare `s et al. 30 Puigjaner et al. 33 S1 Single roll S2 Single roll oriented diagonally S3 Single roll elongated toward two opposite horizontal edges S4 Nearly toroidal roll S5 Four roll structures, each one with its axis perpendicular to one sidewall S6 Two parallel rolls S7 Structure S3 with merged ascending currents S8Two asymmetric counter-rotating rolls aligned along one of the x ¼ z diagonals (similar to S2) S10Two asymmetric counter-rotating rolls aligned along one of the x ¼ z diagonals, similar to S8 possessing no symmetry element

FIG. 10 .
FIG. 10.Velocity streamlines for Ra x ¼ 8.34 Â 10 4 where (a) the front and back walls of the cavity are set to adiabatic while the other two sidewalls are set to be perfectly conducting, and (b) all sidewalls are perfectly conducting, where the blue color represents a lower velocity and the red color represents a higher velocity at a given point in time.The top part of the panels (a) and (b) correspond to the resting configuration of the flow before and after the convective burst, where the velocity of the fluid is close to zero (i.e., a quasi-stationary state is attained, as explained in the text).

FIG. 11 .
FIG. 11.Graphical depiction of the convective state for the case where the two sidewalls are conducting and the front and back wall are adiabatic, in keeping with the description method of Pallare `s et al.31

FIG. 12 .
FIG. 12. Snapshots of flow streamlines (top) and velocity components across the mid-plane of the cavity (bottom) for the case where the sidewalls are conducting and the front and back of the cavity are adiabatic for Ra x ¼ 8.34 Â 10 4 at t ¼ 0.3P.

FIG. 14 .FIG. 15 .FIG. 18 .
FIG. 14. Top view of the cavity showing the two centerlines of the zx plane (dashed) over which the velocity profiles are taken.

FIG. 20 .
FIG. 20.Streamlines colored by velocity magnitude, across one period of vibration for Ra x ¼ 8.34 Â 10 5 and the case where all walls are adiabatic.
28(b), 28(c), 29(e), and 29(f)].It is noticeable, however, that the horizontal velocity U x sampled over the z axis [Figs.28(e) and 28(f)] and the horizontal velocity U z sampled over the x axis [Figs.29(b) and 29(c)] show a velocity which is close to 0, indicating that a possible symmetry axis as shown in Figs.16(a), 16(b), 19(a), and 19(b)

FIG. 21 .
FIG. 21.Streamlines colored by velocity magnitude, taken during the sub-period when the flow embodies a quasi-stationary state for Ra x ¼ 8.34 Â 10 5 , and (a) the case where all walls are adiabatic, (b) the case where two sidewalls are conducting and the other two sidewalls are adiabatic, and (c) the case where all sidewalls are conducting (the blue color represents a lower velocity and the red color represents a higher velocity).

TABLE I .
Comparison of current results with those obtained with the in-house code for the fully conducting 2D case Ra x ¼ 4.00 Â 10 4 , X ¼ 5, Pr ¼ 7, grid size 100 Â 100 (Numean indicates the time-averaged value of the Nusselt number).

TABLE II .
Fluid parameters.