Two-dimensional vibrationally driven solid particle structures in non-uniformly heated fluid containers

For years, the problem of particle accumulation in fluid systems has been of great interest to engineers and scientists due to its ubiquity in natural phenomena and industrial applications.^1–14 In this work, attention is paid to the convective transport of solid particles in a fluid due to pure thermovibrational convection. For the convenience of the reader who may not be an expert in these fields, in this introduction, first, the characteristics and main properties of this unique type of flow are briefly reviewed and then, a short description of the general factors that influence the formation of particle structures in liquids is provided. Finally, recent studies treating particle accumulation driven by thermovibrational effects are discussed, together with a critical assessment of the related gaps to “be filled” (which constitute the natural pre-requisite for the definition of the objectives of the present study).

The simplest way to undertake such a treatment is to consider that thermovibrational convection may be regarded as a variant of natural convection induced by the simultaneous application of a temperature gradient and vibrations. A temperature inhomogeneity in a fluid can give rise to a change in its density, which, together with an acceleration field, typically results in the excitation of buoyancy forces and the generation of fluid motion. In particular, on Earth, the two classic scenarios leading to the onset of buoyant natural convection can be articulated into two main variants, namely, (a) the case where the temperature gradient is perpendicular to the driving force (leading to the so-called Hadley flow¹⁵) and (b) circumstances where the temperature gradient is parallel to the driving force (a very classical problem in fluid-dynamics, generally known as Rayleigh–Bénard or RB convection¹⁶). In standard terrestrial phenomena or technological applications, gravity constitutes the driving force leading to these two modes of convection; however, if it is replaced with (or superimposed onto) vibrations, thermovibrational flow is enabled. Seminal studies pioneering this type of fluid motion trace back to the 1970s and 1980s when Gershuni and co-workers^17–19 first investigated the effects of vibrations on a fluid in weightlessness conditions.

Before delving into a brief review of the physics of this type of flow, it is also worth mentioning that scientist's interest in these phenomena coincides not by chance with the advent of microgravity platforms and the possibility to conduct experiments in space. More specifically, the influence of vibrations on fluid systems in space has evolved over the years from a curiosity driven problem to one attracting significant attention due to its connection with the undesired effects that the appreciable levels of residual acceleration present on space fairing vehicles and orbiting platforms, otherwise known as g-jitters, can have on fluid systems.^20–23 These disturbances are caused by aerodynamic and aeromechanical forces and take the form of high frequency, low amplitude time-periodic displacements. The interested reader may consult Nelson²⁴ for a detailed characterization of these effects in relation to material processing in space.

Given these premises, thermovibrational flow should be treated as a standalone type of convection. Leaving aside for a while the related technological implications, the complexity of this specific problem is intriguing to scientists as many factors contribute to its solution(s).^17–19 These include, among others, the type of fluid considered (Newtonian/non-Newtonian, see, e.g., Refs. 25–27), the direction of the vibrations with respect to the temperature gradient,^28,29 the boundary conditions of the system (including the heating conditions at the walls³⁰ and the shape and size of the fluid container³¹), and the shaking amplitude and frequency. Accounts of the specific effects of the individual and combined factors mentioned above can be found in existing books, see, e.g., Ref. 32, and are not examined in detail in this work so as not to increase excessively its overall length.

When particles are introduced into such convective systems, their response to the fluid flow follows standard rules already known for other fluid–particle systems; namely, their behavior can depend on inertial mechanisms and other factors causing the particles to react in a more or less timely way to external stimuli. Delays in such a response can induce a departure of particle trajectories from the streamlines of the carrier flow, thereby causing non-linear effects that accumulate in time. These processes have been investigated to a large extent for a variety of flows resulting in a set of well-established results and theories.^33–39 

When thermovibrational flows are considered, however, additional unique properties about the interaction between the fluid and particles stem from the “cyclic” nature of the considered flow itself. Thermovibrational flow generally oscillates in time with a frequency that is, in most cases, proportional to that of the imposed shaking (identical, a multiple, or a sub-multiple), and this allows the aforementioned non-linear fluid–particle interplay to undergo an iterative or self-sustained process, which over recent years has been found to result in previously unknown phenomena. Indeed, particles have been observed to cluster and form very regular structures with complex two-dimensional or three-dimensional morphology,^40–46 which resemble that of geometrical entities, such as the “conics” and “quadrics,” whose properties are the main subject of analysis in the field of projective geometry.

While the hidden mathematical laws underlying the emergence of such “order” have not been identified yet, much emphasis has been put on determining the critical links between the morphology of these aggregates (i.e., their shape as seen by an external observer) and the aforementioned “environmental factors,” i.e., the shape of the container hosting the fluid–particle mixture, the directions of vibrations with respect to its walls, and the thermal behavior of such boundaries.^41–44 All of these studies have assumed walls at a constant uniform temperature in order to support a unidirectional thermal gradient (needed to produce the thermovibrational flow itself).

Some effort has also been provided in an attempt to classify these phenomena in terms of spatial multiplicity of the underlying attractors, i.e., the number of loci where the particles tend to be collected with time. These have been found to manifest themselves always in “couples”; i.e., although in some circumstances, coalescence of the emerging structures is also possible and these may display a strong spatial asymmetry, particles always accumulate in two distinct regions of the physical volume inside the container.

In the present work, this specific problem is further explored by replacing the simple thermal boundary conditions assumed in earlier studies with more complex ones. The main objective is a more in-depth analysis of the multiplicity problem when the aforementioned main simplifications, i.e., the unidirectional nature of the imposed thermal gradient, are removed.

By doing so, we pursue a more complete characterization of this recently discovered class of phenomena. Leaving aside for a while the purely theoretical implications that attach to such a purpose, the present work also aims to discern additional useful details, which in the future may lead to the definition of new vibration-based “control” techniques for the contactless manipulation of particles dispersed in a fluid. The related technological implications would be of great value. Suffice it to say that the properties of many composite/multiphase materials depend on the effective distribution of the related minority phase into a hosting matrix (the majority phase). Having the ability to control these would lead to significant improvements in many fields.^9–14 

The problem is tackled under the constraint of two-dimensionality (2D flow). Such approximation is chosen primarily for numerical convenience. We wish to remark that in earlier works on this specific subject (see, e.g., Ref. 41) where three-dimensional (3D) simulations were carried out for cubic cavities with thermal boundary conditions being kept “unchanged” along the spanwise direction, the emerging particle structures were found to be the “projection” in 3D of those found in 2D (i.e., cylindrical surfaces with axes parallel to the spanwise direction and cross-sectional shapes identical to that obtained in the framework of 2D studies). As three-dimensional effects have been reported to be mild in those cases (a departure from the morphology found in 2D becoming evident only in proximity to the walls delimiting the container in a direction perpendicular to the spanwise direction), here, 2D computations are intentionally used as a workhorse to investigate if and how the multiplicity of particle attractee revealed by earlier analyses can be somehow altered by removing the unidirectional temperature gradient constraint (which makes the two-dimensional approximation a logical choice in the context of the present investigation).

The vibrations are modeled mathematically as a sinusoidal displacement in time, characterized by an amplitude b (m) and an angular frequency ω (rad/s), where ω = 2πf 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) can then be formally obtained by taking the second derivative of Eq. (1),

where

Here, two new cases are considered, namely, (1) the situation where vibrations are perpendicular to the adiabatic walls of the square container and the temperature gradient effective through the fluid can change its sense through the cavity due to a non-uniform temperature distribution along the walls parallel to vibrations [case (1)] and (2) a configuration where vibrations have an oblique direction and the corners of the cavity are thermally controlled in such a way that they have different temperatures if pertaining to the same wall [opposed corners having the same temperature, case (2)].

As evident in Fig. 1, by assuming a reference system with origin in the left bottom corner of the cavity, the vibrations are set along the x axis in case (1) and along the direction x = 1 − y in case (2).

The balance equations for mass, momentum, and energy under the assumption of incompressible flow read

where $V_=[u,v]$ and p and T are the problem primitive variables, i.e., velocity, pressure, and temperature, respectively. Moreover, ρ is the fluid density, μ is the fluid dynamic viscosity, and α is the fluid thermal diffusivity. Using the Boussinesq approximation, the density can be considered everywhere constant with the exception of the momentum production (buoyancy) term at the right-hand side of Eq. (5), which can be expanded as

where T₀ is a reference temperature, ρ₀ is the density for T = T₀, and β_T is the so-called thermal expansion coefficient by which variations in the fluid density are directly linked (via a linear relationship) to the gradients of the temperature in the fluid. By substituting this term into the momentum equation, assuming that the individual term $ρ0g$ is absorbed into the modified pressure term and scaling all lengths by the size of the cavity (L), the velocity by (α/L), the time by (L²/α), and the pressure by (ρ₀α²/L²), all these equations can be cast in a 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., ν = μ/ρ₀). Moreover, Ra_ω appearing in the buoyancy term is the vibrational Rayleigh number, analog to the classical Rayleigh number used in standard gravitational convection problems (Ra), namely,

Last, Ω is the non-dimensional angular frequency of the vibrations defined as

where P is the nondimensional period of vibrations.

The initial conditions for each case considered in the present work read

Figure 1(a):

Figure 1(b):

In the light of the information provided in Sec. II A, the boundary conditions for the solid walls read

Case 1 [Fig. 1(a)]:

Case 2 [Fig. 1(b)]:

With the introduction of suspended particles into the liquid comes the need to appropriately take into account the multiphase nature of the considered problem. Here, we adopt an Euler–Lagrangian framework, namely, fixed Eulerian grid modeling of the carrier flow, with the motion of suspended particles being detached and tracked by means of a Lagrangian approach (where, however, the forces exerted by the fluid on the particles are adequately accounted for, a numerical strategy generally known as “one-way coupling”).

Accounting for the associated particle properties (density, shape, and size), obviously, leads to the definition of some additional non-dimensional groups. These are the density ratio (of particles to fluid) denoted by ξ, along with the particle Stokes number arranged into a form depending on the particle radius and the characteristic length of the cavity,

Moreover, it is implicitly assumed that all particles within the system are identical and perfectly spherical in shape. Accordingly, the Lagrangian (Maxey–Riley) equation for particle tracking can be cast in a compact form as

where the additional non-dimensional parameter γ reads

and it accounts for the non-dimensional amplitude of the acceleration induced by vibrations. Moreover, V_p = [u_p, v_p] is the particle velocity and Re_p is the related instantaneous Reynolds number, defined as

where f(Re_p) is a corrective factor required to account for the departure of the drag from the classical Stokes law (Clift et al.⁴⁷),

As explained previously, with this model, the fluid within the considered square container affects the motion of particles, but its dynamics are not altered by the presence of particles. In general, this can be considered a reliable (realistic) assumption if the considered fluid system is “dilute,” i.e., the space between particles is considerably greater than their diameter. Still relying on this hypothesis, particles within the domain are assumed to never interfere with each other. This implies that the total number of particles visible in the figures presented in this study is used solely for the purpose of aiding the visualization of the particle “attractors,” which in the framework of a one-coupling strategy behave as “undisturbed” templates for the accumulation of the dispersed solid mass. The interested reader is referred to the recent work by Lappa⁴⁶ for a more exhaustive elaboration of these interesting concepts and the ranges of applicability of the various multi-way coupling approaches existing in the literature to the specific case of thermovibrationally driven particle accumulation structures. Here, we limit ourselves to mentioning that, from a mathematical/modeling standpoint, it is also essential that the value of Stokes number is ≪1 in order to meet the conditions for which the aforementioned Maxey–Riley equation is valid.

Initially, particles are distributed consistently throughout the cavity, and in the event where a particle approaches a wall, it is forced to maintain a distance from it not smaller than its radius (accordingly, it will sit or slide along the boundary until its directional velocity component perpendicular to the boundary changes or the vibrationally induced acceleration forces it to leave the wall).

The numerical simulations included in this work have been produced by solving the balance equations presented in Eqs. (4)–(6) in the framework of a finite volume method for incompressible flow. In particular, the computational platform ANSYS Fluent has been used. This software relies on the well-known PISO algorithm (Pressure Implicit Split Operator) to ensure adequate pressure–velocity coupling. The details of this time marching procedure, not reported here for the sake of brevity, can be found in a number of works previously published by the present authors (see, e.g., Refs. 29, 30, and 46). It relies on a precise hierarchy of computational steps, where, first, a simplified version of the momentum equation (obtained by dropping the pressure term) is integrated in time to obtain an “intermediate” velocity field and then a Poisson equation is solved for the determination of the previously neglected pressure. This Poisson equation formally stems from expressing the physical velocity as the summation of the intermediate velocity and the pressure gradient and forcing this linear combination into the continuity equation. This system of equations is typically complemented by boundary conditions represented by the effective velocity on the boundary of the physical domain; accordingly, Neumann boundary conditions are used for the pressure. With ANSYS Fluent, the PISO algorithm is implemented in the frame of an implicit temporal approach. For the present simulations, a second order upwind spatial discretization scheme has been adopted for the convective terms appearing in the momentum and energy equations (standard central differences have been employed for the diffusive terms only).

For what concerns the integration of the particle transport equation [Eq. (20)], ANSYS Fluent relies a 5th order Runge–Kutta scheme. Unlike the approach implemented for the fluid phase, in this case, only the particle velocity appearing in the drag contribution is dealt with in an implicit way from a temporal standpoint, whereas all the other terms are treated explicitly. For the sake of completeness, we should also mention that, in order to couple the Eulerian and Lagrangian kernels, we have used the “barycentric method” to determine the fluid velocity at the instantaneous positions occupied by particles, and a bi-linear interpolation and first-order Euler schemes for the material derivative DV/Dt [both V and DV/Dt, indeed, appear at the right-hand side of Eq. (20) and need to be determined from the corresponding Eulerian fields accordingly].

Given the unavailability of results by other authors about thermovibrationally driven particle–fluid mixtures, following the same approach already undertaken by Ref. 46, we have assessed the accuracy of the numerical solver described in Sec. III through consideration of an in-house benchmark problem; i.e., the results provided by two independent computational platforms have been compared.

For this two-dimensional test case, vibrations perpendicular to the direction of the temperature gradient (Hadley-flow-like configuration) and adiabatic conditions at the sidewalls have been considered. The related non-dimensional parameters are summarized in Table I. A mesh size of 90 × 90 has been used.

The agreement found between the solver originally developed by Lappa⁴¹ and the present one (ANSYS Fluent version 2020R2) has been verified by looking at both the x and y velocity components, signals being taken at three locations in the cavity P1, P2, and P3 where P1 [0.25, 0.75], P2 [0.25, 0.5], and P3 [0.25, 0.25].

As quantitatively substantiated by Fig. 2, almost no distinction can be made between the results produced with ANSYS Fluent compared to those obtained with the in-house code by Ref. 41 (which unlike ANSYS Fluent is completely based on an explicit-in-time approach). An additional level of validation obviously stems from comparison of the emerging particle structures in terms of patterning behavior and size (Fig. 3). The minor differences in terms of lobe size visible in Fig. 3 can be ascribed to the different interpolation schemes used by Ref. 41 and ANSYS Fluent (the reader being also referred to the arguments reported in Ref. 46). We conclude that the “physics” underlying the considered phenomena is captured with a similar level of fidelity by both codes.

The validation study presented in Sec. IV is naturally complemented here by a mesh refinement assessment conducted to verify that in addition to the “physics,” also, numerical aspects have been finely tuned. Just like in classic thermal convection problems, the parameter known to influence the “strength” of convection is the Rayleigh number [shown in Eq. (12)]; in thermovibrational problems, this role is played by the vibrational Rayleigh number. The higher the Rayleigh number, the more complex the flow response and, therefore, the greater the number of grid points required to accurately capture the underlying physics. As in this study, the vibrational Rayleigh number is fixed to Ra_ω = 1.56 × 10⁴ and Pr = 8; accordingly, the mesh refinement study is carried out assuming these specific values (see Tables II and III).

The outcomes of the grid refinement analysis are shown in Table IV. The x and y components of the velocity refer to a location P where the velocity magnitude of the flow attains its highest value (P [0.035, 0.08]).

No observable differences occur upon an increase in the mesh resolution from 90 × 90 to 105 × 105; therefore, the 90 × 90 mesh is selected.

Following a logical approach, where the simplicity of the thermal boundary conditions considered in preceding studies is gradually replaced by more involved thermal configurations, first, we examine case (1) [Fig. 1(a)]; i.e., the vibrations are imposed along the x-direction and three alternating temperature gradients are set in the y-direction. The non-dimensional groups (γ) and (St) are varied parametrically over a range of 1.5 × 10⁶–6.5 × 10⁷ and 1.39 × 10⁻⁵–8 × 10⁻⁴, respectively. Since earlier efforts^31,41 have shown that particle accumulation does not occur in the limit as the density of particles becomes equal to that of the surrounding liquid, here, a value of ξ > 1 is considered (without loss of generality, we fix it to ξ = 1.85, the reader being referred to the cited literature for a detailed information about the influence of the density ratio on the considered phenomena⁴⁴). Since in the present conditions, many structure types are observed, most conveniently, the following nomenclature is used to adequately characterize the new findings; i.e., the particle accumulations are referred to as the “four-loop” type, the “central pillar accumulation,” the “two central structures,” and the “two side extensions.”

In this subsection, a brief description of each of the abovementioned four realizations is elaborated starting from the “four-loop” type structure, identified by the circle symbol (⁠$∙$⁠) in Fig. 4. This structure type has been found at three locations in the map, i.e., at medium to high values of both γ and St. Moreover, as a distinguishing mark with respect to other cases, symmetry holds with respect to the x axis for y = 0.5 and about the y axis for a time dependent value x_t, where x_t is the point along the x axis representing the location of the vertical centerline of the structure. Four “loops” (or particle-dense “circuits”) can clearly be recognized, each resembling a Cartesian oval (egg shape). From a practical (interpretative) standpoint, these ovals may be regarded as “attractors” to which all particle trajectories tend, easily revealed by numerical simulations when a statistically relevant number of particles is considered. In this regard, we wish also to recall that, as originally pointed out by Lappa,⁴¹ after a transient time, the related particle formations attain a well-defined (asymptotic) shape and size. Although the instantaneous structures keep oscillating along the direction of the vibrations, they give the observer the illusion of a “solid body” moving back and forth (which means that the instantaneous structures reported in Fig. 4 are fully representative of the related time-averaged counterparts).

The next formation we lend our attention to is the “central pillar,” found at six instances, represented by the square symbol (■) in Fig. 4. As qualitatively substantiated by the insets in this figure, only one bulk structure emerges in this case, in contrast to the other three types where distinct particle aggregates exist, as also made evident by the “interposed” clear (particle-free) fluid. The central pillar is formed primarily due to wall effects. Indeed, the simplest way to justify its existence is to consider that the displacement of the cavity along the x axis (i.e., the shaking direction) causes the particles to cluster along the vertical sidewall and detach when the cavity is displaced in the opposite direction. These behaviors have been found to be a key factor in influencing particle accumulation in thermovibrational problems when considering differentially heated cavities (Hadley flow type configuration) in both laminar⁴¹ and turbulent conditions.⁴⁵ The interaction with the walls, however, is not the only influential factor. Convection still plays a role as witnessed by the tendency of the structure to contract in an hourglass shape as time passes (where the smallest extension along x is visible at y = 0.5).

The next case of the sequence, i.e., the accumulation represented by the diamond symbol (⁠$⧫$⁠) in Fig. 4, is the rarest of all formations. Indeed, only two instances have been detected (for the range of parameters considered here). Denoted as “two central structures,” this variant shows up in the region of medium values of γ and high values of St. This case bares much similarity to the (■) morphology; indeed, an even more pronounced constriction is visible at y = 0.5. What sets it apart from the other realizations is that the central column splits into two structures and a finite layer of clear fluid emerges through the center of the overall particle-dense region. The perimeter of the structure is also more defined when compared to the central pillar accumulation (■), where the structure is “closed” near the top and bottom walls.

Finally, we turn to the “two side extension” variant, represented by a triangle symbol (⁠$▴$⁠) in Fig. 4. The central region of accumulation visible here is very similar to the aforementioned central pillar represented by the square symbol (■). In this case, however, the pillar is ornated by two elongated ovals at the left and right sides.

Having completed a description of the patterning behavior, as a concluding observation for this subsection, we wish to remark that two global symmetrical trends apply to all these cases. Indeed, symmetry can be observed about the x and y axis midplanes for all the formations.

This subsection is devoted to a quantitative assessment of the sensitivity of these structures to the problem parameters in their respective ranges of existence. In this regard, we start from the remark that while a decrease in the particle size (i.e., a decrease in St) is expected to lead to less compact structures, a higher vibrational amplitude γ should have the opposite effect.⁴⁴ It has also been demonstrated that under certain conditions, the tendency toward more cohesive formations due to an increase in the value of γ may be offset by a decrease in St and vice versa.

Notably, a similar scenario can be discerned in Fig. 4, where diagonal trends in aggregate types are evident from the N-W to S-E regions of the map when γ decreases and St becomes larger. These observations are naturally complemented by the findings collected in Fig. 5, where the “characteristic length” of the structure has been reported as a function of time for each particle accumulation type. Given the different morphologies and the intrinsic multiplicity of the different formation types, the characteristic length is defined as follows: “four-loop” type [Fig. 5(a)], the maximum size of a single loop along the y direction; “two central structures” [Fig. 5(b)], the maximum size of one of the two pillars along the x direction; “central pillar accumulation” [Fig. 5(c)], the horizontal width of the “necked” region at half height of the cavity; and “two side extensions” [Fig. 5(d)], the maximum size of the side loop along the y direction. In particular, in each panel of Fig. 5, the evolution of two cases, differing in terms of γ and St, has been reported. The most interesting outcome of this figure resides in its ability to show that, in general, by compensating a decrease in St with an increase in γ, the final size of the structures (and the related formation time) can be kept unchanged. By contrast, as qualitatively and quantitatively substantiated in Fig. 4, significant changes (simultaneous increase or decrease) in both St and γ can cause the transition from one structure type to another.

Another influential parameter is considered here (already identified in the frame of earlier studies where simpler thermal boundary conditions were considered), i.e., the frequency of vibrations (Ω) defined in Eq. (13). This parameter was previously fixed to Ω = 2 × 10³. In this section, three additional (lower) values of Ω are considered for the following structure types: “two side extensions” (⁠$▴$⁠, γ = 1.5 × 10⁷ and St = 8.01 × 10⁻⁵) and the “four-loop” type structure (⁠$∙$⁠, γ = 4 × 10⁷ and St = 1.39 × 10⁻⁴).

The corresponding modifications in the “two side extension” variant (for Ω = 1.2 × 10³, Ω = 1.3 × 10³, and Ω = 1.5 × 10³) are summarized in Fig. 6. Interestingly, with an increase in the frequency of vibrations, the side ovals become more extended along the y axis. The central pillar, however, remains undisturbed up until a value of Ω = 2 × 10³, where two main lobes emerge resembling those pertaining to the “two central structure” case. This reinforces the idea of a seamless transition between different structure types as already discussed to a certain extent in Sec. VI A 2.

The quantitative increase in extension along both the x and y axis is detailed in Figs. 7(a) and 7(b). The structure formation time (defined here as the time required to obtain a state where the characteristic size of the structures and their distances remain constant) is also included in Fig. 7(c) for the sake of completeness. A rise of the formation time is observed when Ω is increased, and these findings align with those by Ref. 44.

The analogous results for the “four-loop” type structure are depicted in Figs. 8 and 9 (the range of values investigated here includes three additional values of Ω: Ω = 1.6 × 10³, Ω = 1.7 × 10³, and Ω = 1.85 × 10³).

A similar trend can be recognized in this case with a global increase in both the length of the x and y extensions of the ovals, together with a growth of the formational time when Ω becomes higher.

All these findings can be naturally combined with those pertaining to the 2D map shown in Fig. 4. The outcomes of this exercise, by which the original 2D map is expanded through consideration of an additional degree of freedom (i.e., the vibration angular frequency), are shown in Fig. 10.

In this section, the second configuration depicted in Fig. 1(b) is considered [case (2)]. As already explained in Sec. II A, the N-W and S-E corners of the cavity are set to T_cold and the N-E and S-W corners to T_hot, leading to a symmetry about the x = y line.

Before starting to deal with this situation, at the risk of repetition (see Sec. I), we wish to recall once again that the direction of vibration is known to be an important influential factor in relation to the response of thermovibrational systems.^41,42 Here, in particular, for consistency with the results presented for the companion case with non-inclined vibrations, the following case is assumed as a “basic” configuration (to assess the changes induced by a variation in the shaking direction): Ra_ω = 1.56 × 10⁴, Ω = 2 × 10³, γ = 6.5 × 10⁶, St = 3.04 × 10⁻⁴, ξ = 1.85, and ϕ = 0.

As shown in Fig. 11, the emerging structure bares much resemblance to the “two side extension” variant [the triangle symbol (⁠$▴$⁠) in Fig. 4]. However, some minor differences can be spotted. First, the central pillar does not dispose of a contracted mid-section and appears to be of a consistent width from 0 < y < 1. Second, if one looks at the location from which the side loops originate, while for the case presented in Sec. VI A, the ovals originate from a stream of particles attached to the bottom wall (y = 0, where the wall segments were set to T_hot), for the corner-heated cavity, these originate from the heated wall sections.

The zoo of solutions produced when the inclination angle is varied is finally summarized in Fig. 12 (over the range 0 ≤ ϕ ≤ 3π/4). It can be seen that when ϕ is changed from 0 to π/2, the entire pattern simply undergoes a complete 90° rotation. For a relatively small departure from these limiting conditions (e.g., ϕ = π/12 or π/6 and 5π/12 or π/3, respectively), the central column begins to warp in the direction of the hot corners.

Finally, a range of values of ϕ exists (2π/9 ≤ ϕ ≤ 3π/10) where four distinct loops are formed. Owing to the symmetries embedded in the problem, the scenario for ϕ spanning other sub-ranges of the 0 ≤ ϕ ≤ 2π interval can be inferred on the basis of similar arguments (see, e.g., the inset for ϕ = 3π/4).

Placing the present results in a wider theoretical context requires a short excursus on the unique properties of this class of particle attractors and what has been understood until now in terms of underlying cause-and-effect relationships. In this regard, as already discussed to a certain extent in Sec. I, it is worth recalling that the present work adds another piece to the puzzle related to the existence of particle attractors in vibrationally driven flows, which do not depend on particle interactions (whose first piece was placed by Lappa⁴¹ with the discovery of a new class of self-induced particle aggregates in vibrated non-isothermal fluid systems in the framework of a one-way coupling approach). From looking back at these previous works, some useful arguments can be drawn, which can help to interpret the present one and elaborate a more general “view” of these phenomena and useful generalizations.

In particular, compelling evidence has been provided^31,41 that particle attractors are made possible by the interplay of inertial and thermovibrational effects as indirectly proven by the suppression of particle structures in the limit as the density of the particles becomes equal to that of the surrounding liquid or the thermovibrational effect is disabled (while still allowing non-iso-dense particles to feel the acceleration produced by vibrations^31,41). Convincing arguments have also been elaborated about the direct relationship between the effective spatial shape of the “attractors” (directly connected to the morphology of the particle structures emerging in the framework of one-way coupled numerical studies) and the relative direction of vibrations with respect to the solid boundary of the systems.⁴¹ According to these studies, if the angle between vibrations and the walls of the system is changed while vibrations remain perpendicular to the imposed temperature gradient, a zoo of structures can be obtained. Subsequent efforts⁴² have revealed that further modifications in the morphology can be induced by allowing vibrations to change their inclination with respect to the imposed temperature gradient, which has led to the main conclusion that the effective physical realization in space of particle structures is dictated by the threefold relationship between the shaking direction, that of the imposed (unidirectional) temperature gradient and the effective shape of the container physically hosting the fluid–particle mixture.

Other works of relevance to this subject have provided evidence that additional levels of complexity can be brought in by symmetry breaking effects emerging in the carrier flow itself in a “spontaneous” way (i.e., due to an intrinsic bifurcation of this flow, which would occur even in the absence of particles⁴⁴) or produced as a result of the back influence of dispersed solid matter on such a flow (i.e., owing to a two-way coupled mechanism⁴⁶).

Given these premises, the major contribution of the present study should, therefore, be seen in the new findings about the multiplicity of the attractors, i.e., the number N of attracting loci coexisting in the physical space at the same time. Evidence has been provided that this property can be somehow controlled by acting on the degree of thermal inhomogeneity of the fluid.

Through a generalization of a large number of results, it was previously concluded that this multiplicity N could not exceed a value N = 2 as the particle structures were always manifesting in couples or as a single large cohesive unit produced by the coalescence of two initially distinct aggregates. Although conducted under the constraint of a 2D flow, the present simulations have shown that even if the analysis is still limited to a case as simple as a square cavity, conditions can be identified for which several distinct formations (with N > 2) can co-exist within a single fluid domain.

In particular, by indicating with M the number of inversions in the temperature gradient through the fluid, the multiplicity for the first setting considered in the present work [case (1) shown in Fig. 1(a)] can be expressed as follows:

The significance of this inequality, which holds for a fixed value of ξ, Ra_ω, and Ω (1.85, 1.56 × 10⁴, and Ω = 2 × 10³, respectively, in the present work), can be further explained or elucidated as follows: by changing the amplitude of vibrations (γ) for a fixed size of particles (St) or, vice versa, changing the Stokes number for a fixed γ, topological changes can be induced in the emerging particle pattern, with the number of possible realization constrained to be less than or equal to 4.

For the second paradigm considered, i.e., the cavity shown in Fig. 1(b) [case (2) for which M = 1 and the vibrations can form an angle ϕ with the x axis], this relationship still holds provided a corrective factor is introduced to account for the inclination of the shaking direction, i.e.,

which indicates that the relative direction of vibrations can also affect the multiplicity parameter.

The problem related to the formation of particle accumulation structures in thermovibrationally driven systems still carries a number of interesting questions, some of a general nature, other more system specific. In the present work, a first attempt has been made to understand if and how behaviors and mechanisms known to be effective under unidirectional temperature gradients can be mapped into systems with non-uniformly heated boundaries.

Emerging structures have been categorized into specific classes depending on their appearance in the physical space. After running a statistically representative number of numerical simulations, it has been demonstrated that a colorful spectrum of variants is possible, which differ in terms of topology, morphology, and number of attracting loci. For a fixed couple (ξ, Ra_ω), these accumulations exist through a multi-dimensional space of parameters [particle size (St), amplitude (γ), frequency (Ω), and angle (ϕ) of vibrations] as a result of the reverberation of physical processes driven by thermovibrational (convective) and inertial (particle-related) effects into a confined space. The major outcome of the present study is that such phenomena can also be “tuned” or “controlled” to a certain extent by acting on the degree of thermal inhomogeneity of the system. On increasing the number of inversions of the temperature gradient inside the fluid, the multiplicity of the emerging structures grows accordingly and transition among different realizations can be obtained on varying γ and/or St. Some correlations or relationships have been introduced in an attempt to model these behaviors. Although these may be regarded as still incomplete realizations of a more complete theory that shall be formulated to predict the properties of particle formations in all situations and as a function of all the influential parameters, this work should be considered a further step toward the complete characterization of the related dynamics. An interesting prospect for the future is the extension of this line of inquiry to the more general three-dimensional case for which the multiplicity of structures and the related response to changes in the system parameters and boundary conditions are expected to become even more intricate.

This work was supported by the United Kingdom Space Agency (STFC Grant Nos. ST/S006354/1, ST/V005588/1, ST/W002256/1, and ST/W007185/1) in the framework of the PARTICLE VIBRATION (T-PAOLA) project.

The authors have no conflicts to disclose.

Georgie Crewdson: Data curation (equal); Investigation (equal); Project administration (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Matthew Evans: Data curation (equal); Investigation (equal); Software (equal); Visualization (equal). Marcello Lappa: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are openly available in Pure at https://doi.org/10.15129/cdd1d3fd-c533-49b4-9197-a24b05bc1e3e.