Characterization of two-way coupled thermovibrationally driven particle attractee

Fluids with dispersed particles play an important role in numerous contexts. They really stand at the crossroad of fundamental physics,^1–3 astrophysics,^4,5 chemistry,⁶ biology,^7,8 materials science,⁹ and thermal and mechanical engineering.^10–12 Related problems include (but are not limited to) fluid-dynamic instabilities in complex fluids,^13–15 the interaction of an immiscible phase with a turbulent flow,^16–19 particle ordering in laminar convection,^20–33 and the fascinating prospect of using some of these phenomena to produce “active” matter, i.e., parts or components that can recognize and bind to each other and form well-defined networks or systems.^5,9 Other intriguing directions of research originate from the possibility to interpret theoretically all these behaviors as they were the manifestation in the physical space of complex non-linear dynamics governed by overarching “attractors” in the phase space (“attractee”). The advantage relating to this point of view obviously lies in the ensuing process of abstraction, by which specific cases (like the one addressed in the present work) can be used as a testbed for the definition of more general problems or fundamental questions.

In such a context, it is worth highlighting that much attention has been attracted by turbulent flows given their ability to promote particle clustering via different possible routes and because of the myriad technological applications they are relevant to. The underlying mechanisms have been studied over many years (see, e.g., Ref. 17 and references therein) together with the tendency of these flows to induce diametrically opposite behaviors, i.e., particle scattering (see, e.g., Ref. 34). Although the variety and dispersion of these results still call for an effort to unify the information to general criteria, more recently, some physicists and engineers have placed their hopes in laminar (deterministic) flows as possible candidates for the identification of new particle attractors. Notably, this specific line of inquiry has succeeded in showing that remarkable exemplars can be found in the category of simply time-periodic flows.^5,24,35,36

As the reader might have realized at this stage, obviously, the final goal of all this research is a more exhaustive (academic) characterization of all these phenomena (and, indeed, relevant theoretical explanations have been elaborated for the ability of particles to self-organize in both turbulent and laminar cases). In particular, it has been clarified that a large part of these mechanisms rely on purely topological relationships. While point tracers should, in general, be considered as mere abstractions, which do not reflect physical reality, real particles have finite size and often different density with respect to the fluid hosting them. For these reasons, their trajectories and the streamlines of the carrier flow do not overlap exactly. Since the small departures caused by particle drag, inertial effects, and other body forces present in the considered system tend to accumulate, the distribution of particles can undergo significant changes as time elapses. This process can finally result in regions of clear fluid (with relatively small or negligible presence of solid matter) and, by contrast, dense clusters formed by the spontaneous accumulation of particles. Notably, as this behavior somehow resembles that of a gas undergoing changes in the local density due to significant gradients of velocity (this being often the case in typical problems of compressible aerodynamics), the evolution in time of finite-size particles can be assimilated to that of a compressible medium (Refs. 24, 32, and 33 and references therein).

A typical fingerprint of this class of events in laminar flows is their extremely ordered appearance (as opposed to the fractal morphology of the disordered aggregates induced by turbulent flow). These patterns can display a variable number of dimensions depending on the considered circumstances. As an example, they are essentially one-dimensional (manifesting as a thread of particles) when the time-periodic nature of the carrier flow is due to the propagation of a wave along the axis of a basic roll,^{21,22,24,37,38} or two-dimensional if two counter-propagating waves are present. In the latter case, particles generally accumulate in those regions of the space where the amplitude of disturbance resulting from the combination of the two waves is negligible.^31,39,40

Much more complex shapes can be produced if the time periodicity of the flow is due to other effects, and this is indeed the case of thermovibrational convection, i.e., fluid motion induced in a non-isothermal liquid by the application of vibrations.^41–55 Interestingly, the study of particles dispersed in a fluid undergoing this type of convection inside a compact cavity has revealed a zoo of particle structures with surprising geometrical perfection (such as cylinders, hyperboloids, ellipsoids, and conical surfaces^5,56–59), which suggests all these dynamics might be governed by hidden (not yet clarified) precise mathematical laws.

Existing numerical works on the subject have added several pieces to this puzzle. It has been understood that the underlying mechanisms are essentially independent of inter-particle forces, i.e., particle-to-particle interactions are not required to justify the existence of attractors. This should be regarded as a notable difference with respect to other aggregates where vibrationally induced hydrodynamic forces play a crucial role (such as those studied, e.g., by Refs. 60–63). Attractive forces are not needed because the existence of accumulation loci, all having the shape of “quadrics,” is supported by a delicate phenomenon of synchronization between the inertial response of the dispersed particles to the variable acceleration produced by vibrations and the thermovibrational convective effect per se.^56–59 Most remarkably, this specific property allows the emergence of particle structures also in relatively dilute dispersions.

As a fleeting glimpse into the existing literature would immediately reveal, the umbrella of ongoing research on this specific subject has essentially drawn on studying the influence of asymmetries present in the considered physical system on this synchronization process (asymmetries due to the shape of container,⁵⁷ due to the direction of vibrations with respect to the walls of the cavity,^58,59 or asymmetries spontaneously produced due to time-averaged convective effects,⁶⁴ or other mechanisms enabled when the carrier thermovibrational flow displays a certain intrinsic degree of turbulence⁶⁵). No analysis has been devoted until now to explore another important question, that is, to what extent these delicate structures can survive in non-dilute fluid–particle mixtures. Given the essentially inertial nature of these phenomena, this is obviously a legitimate question. The volume fraction of dispersed mass should be considered as another relevant influential parameter.

In order to address this problem, the present article presents findings about the back influence that particles can exert on the carrier thermovibrational flow when a relative high “load of mass” is considered (defined as the ratio of the overall solid mass to that of the fluid).

The required theoretical–mathematical two-way coupling framework is elaborated in Sec. II together with a brief presentation of the main objectives of the overall study.

In order to meet the main goal discussed in the introduction, the intended outcomes of the present analysis are set as follows:

• Develop and validate a generalized model for evolving particles in fluid flow described by a hierarchy of mathematical kernels based on the combination of Eulerian (for fluid flow) and Lagrangian (for particles) approaches, by which different degrees of coupling can be implemented for the solid and liquid phases.

• Expand the space of parameters through explicit consideration of the solid mass load (or particle volume fraction) and perform a parametric study aimed at the identification of particle attractors (the “attractee”).

• Determine the regions of this space where, using the attractors as “templates” for the accumulation process, particles can effectively form recognizable structures.

• Evaluate the impact of two-way coupling on the morphology and regularity of these structures (including the possible emergence of chaos or turbulence for some specific conditions).

For consistency with earlier studies on the subject where dilute distributions of particles were examined (under the assumption of one-way coupled flow), here we focus on the same archetypal geometrical configuration considered by Lappa,^56,59 namely, a differentially heated cubic enclosure undergoing vibrations along a direction perpendicular to the imposed temperature difference (as shown in Fig. 1, temperature difference and vibrations are along the z and y axes, respectively, while x plays the role of “spanwise direction”).

By modeling vibrations through a simple sinusoidal function,^47,48 i.e., as $s ¯ lab ( t ) = bsin ( ω t ) n ̂ ¯$ (where $n ̂ ¯$ is the unit vector along the direction of shaking, b is the displacement amplitude and ω = 2πf is the related angular frequency), the time-varying acceleration affecting the entire system can directly be obtained by taking the second time derivative of s(t), that is, $g ¯ ( t ) = g ω sin ( ω t ) n ̂ ¯$⁠, where obviously g_ω is given by the product of b and the square of the angular frequency. These simple considerations lead to the identification of the first two independent problem variables, namely, ω and b. Of course, other degrees of freedom are represented by the nature of the fluid itself and the applied temperature difference.

These initial arguments clarify the concept that, for the phenomena of interest, vibrations can affect particle dynamics via two independent but concurrent mechanisms, i.e., the force which is produced by the time-periodic acceleration field as a results of the different density of the particles with respect to the surrounding fluid (a purely inertial effect) and the flow of thermovibrational nature that tends to transport particles in a more or less efficient way depending on their radius (in turn, as explained before, the latter consists of two contributions, i.e., the aforementioned “mean” and “instantaneous” velocity components). Obviously, another important parameter is the mass load, which (as explained in detail in Sec. II D) depends on the effective number of particles present in the fluid.

A key observation regarding these problems is that experimental measurements aimed to get quantitative data are extremely difficult (almost impossible) in normal gravity conditions, given the disturbing presence of gravity (which would tend to demix particles from the fluid due to sedimentation of flotation). On the other hand, conducting tests in microgravity conditions is extremely expensive. In practice (while waiting for opportunities for space experiments, Lappa et al.⁶⁶), the most convenient way to gain an understanding of the underlying flow physics is to use numerical simulation. This explains why the existing literature on this subject relies to an unusual level on insights obtained via theoretical models and ensuing computations.

In such a context, however, the assumption of one-way coupled particle-fluid systems, which so much success has enjoyed in earlier works on the subject, should obviously be regarded as a non-viable option. For the specific problem considered here, the numerical treatment of these processes shall necessarily imply adequate physical and mathematical modeling of the mechanisms relating to the transfer of momentum from the particles to the surrounding fluid and vice versa. This naturally leads to the aforementioned concept of “two-way” coupling, by which the two considered (solid and liquid) phases can experience a “mutual interference.”

In principle, many approaches could be applied to meet the objectives set at the beginning of Sec. II. These strategies can be split into three main classes according to the nature of the model used to track the dispersed phase: capturing (also known as moving grid or Lagrangian approach, see Ref. 67), Eulerian (also known as fixed-grid approach, see, e.g., Ref. 68) and combined “hybrid” techniques (see, e.g., Ref. 69 for a review) by which special additional markers are tracked separately to reconstruct the interface between the different phases (see, e.g., Ref. 70) or the Eulerian governing equations in the fluid domain are “extended” to the particle domain by constraining the motion inside the particles to a rigid motion.⁷¹ 

Following already successful attempts,^{24,25,28,29,56} here, in particular, we resort to a hybrid formulation by which, while fluid flow is tackled in the framework of a Eulerian description (leading to the classical Navier–Stokes equations valid for a continuum), the transport of particles is handled on the basis of a Lagrangian method that can account for each particle as a separate entity interacting with the “environment.” The novelty of the present work with respect to earlier efforts (Refs. 56–59 and 64) is obviously due to the two-way-coupled nature of the resulting strategy, which (as we will illustrate in Sec. II C) requires proper physical reasoning and modifications or “additional terms” to be considered for the Navier–Stokes equations.

In order to develop further the ideas discussed above, in this section, we lay the foundation of the overall numerical approach by transcending specific cases and introducing the fundamental equations that govern the considered dynamics. Anyhow, for the convenience of the reader, advanced formal discussions are limited to a minimum and the problem is described with a relatively simple formalism (making reference to preceding works when additional layers of knowledge are particularly appropriate).

Moreover, the additional vector quantity S_m at the right-hand side of Eq. (7) is the term required to properly couple the liquid and solid phases (notably, it makes the resulting approach formally similar to that implemented by Ref. 72 to investigate the dynamics of shock waves traveling in a dusty gas, as we shall further detail in Sec. II D).

An important assumption underlying this modus operandi, however, is the condition that the volume of the generic particle is smaller than the volume of the computational cell, i.e., ϕ < 1 where $ϕ = 4 3 π R p 3 / Δ x Δ y Δ z$⁠. In practice, this is formally needed to justify that particles can be treated as microscopic quantities compared with the size of a typical mesh cell. Notably, this requirement may also be regarded as a mathematical constraint for the maximum tolerable number of grid points used to discretize the problem, or alternatively for the maximum possible diameter of the particles for a given grid resolution (the interested reader being referred to Lappa et al.,⁷² for additional elaboration of these concepts).

Other key observations concern the lack of coupling terms to be considered for Eqs. (6) and (8), which also requires proper physical reasoning and justification.

An explanation for the absence of a mass coupling term can be elaborated in its simplest form on the basis of the argument that, while the liquid and the dispersed solid mass exchange momentum (indeed, most of the dynamics presented in this work stem from this effect), there is no corresponding mass exchange between these two phases (the overall mass of liquid and that of particles are conserved separately). To elucidate further the significance of this observation, one may consider that the mass of the present particles is fixed, i.e., they do not undergo any mass decrease or increase in time.

Another key point concerns the thermal expansion coefficient. This is assumed to be one order of magnitude smaller for the solid (e.g., β_T = 2.07 × 10⁻⁴ and 0.255 × 10⁻⁴ K⁻¹, yet at 20 °C, for water and glass, respectively), which explains why the dependence of the solid density on temperature is not taken into account in the present study.

The governing equations must be complemented by the initial and boundary conditions required to close the problem from a mathematical point of view.

Moreover, at t = 0, particles are distributed uniformly and with zero velocity.

For what concerns the particles, in line with earlier studies,^56,57 these are assumed to interact in a non-elastic fashion with walls, i.e., they can approach the solid boundary until at a distance not smaller than their radius is achieved; thereafter, they are allowed to slide along the boundary until the wall-normal velocity component changes sign (allowing the particles to leave the wall).

After building an appropriate mathematical representation of the considered system, the next step obviously relates to the selection of an adequate numerical procedure to solve the underlying equations.

As solution technique we have implemented a projection method based on the use of primitive variables (namely, velocity, pressure, and temperature). The mathematical prerequisite at the basis of this class of methods has its origin in the so-called decomposition theorem for vector fields (the Hodge theorem, see, e.g., Ref. 82). According to this theorem, a general vector field can always be split into a part of given divergence and the gradient of a scalar potential function; thereby, the mathematical problem relating to the need to determine a velocity field that satisfies the continuity equation and develops vorticity (as required by the viscous nature of the fluid) can be articulated in a numerical procedure where these two aspects are formally taken into account separately (in practice, this process is mediated by the relationship between the pressure and the velocity fields as further illustrated below).

As indicated by Eqs. (26)–(28), we have implemented this numerical technique using schemes explicit in time. In particular, given the extremely small velocities involved in the present problem, the convective and diffusive terms appearing in both the momentum and the energy equations have been discretized with standard central differences. In order to improve the coupling between pressure and velocity (see, e.g., Refs. 87 and 88), however, a staggered arrangement has been used for these two quantities, that is, while the components of velocity u, v and w occupy the center of the cell face perpendicular to the x, y, and z axes, respectively (see, e.g., Refs. 89 and 90), pressure is located in the center of each computational cell.

The particle tracking equation re-arranged to have dV_p/dt appearing only at the left-hand side [as made evident by Eq. (17)] has been integrated with an explicit in time 4th order Runge–Kutta scheme.

For all the simulations, the time integration step has been kept smaller than St/Pr in order to resolve properly the particle relaxation time. Moreover, a grid refinement study has been carried out to ensure grid-independence. In line with the mesh assessment conducted by Ref. 50, a resolution corresponding to 80 grid points along each spatial direction has been found to be more than sufficient (this finding reflects the known fortunate absence of boundary layers in thermovibrational flow for the range of values of the Gershuni number considered in the present work).

Given the nature of the problem addressed in the present study, finding relevant results produced by other authors to be used for validation purposes was not straightforward as one would imagine. For this reason, after unsuccessfully sifting through plenty of existing studies available in the literature with different (nonrelevant) foci, we decided to validate the overall theoretical architecture illustrated in Secs. II A–II F and III A through a comparison with the outputs of simulations conducted using a different computational platform (for the same circumstances considered in the present work). As independent software, in particular, we have considered ANSYS Fluent given the availability of a reliable two-way coupled approach for particle tracking directly coded into this package.

For incompressible flow, ANSYS Fluent is based on a primitive-variable method pertaining to the same category of projection techniques described in Sec. III A, namely, the so-called SIMPLE (semi-implicit method for pressure linked equations) strategy (see Ref. 91). In place of the staggered arrangement of variables (described in Sec. III A), however, this software is dependent on a collocated grid approach, which means that all the physical quantities of interest are defined in the center of the grid cells. For this reason (in order to resolve properly the physical coupling of velocity and pressure), ANSYS Fluent relies on the special interpolation schemes by Ref. 92.

We have used this package to solve the classical set of balance equations for mass, momentum and energy reported in Sec. II C using 2nd order upwind schemes for all convective terms (standard central differences have been employed for the diffusive terms only). With ANSYS Fluent, all these terms are typically handled implicitly (with the sole exception of the buoyancy contribution, the only entity being treated explicitly). Moreover, in order to accelerate convergence, this software takes advantage of a classical algebraic multigrid scheme (AMG) with standard parameters (the so-called Gauss–Seidel smoother, see, e.g., Ref. 93).

A dedicated discussion is also needed for the particle motion equation. Fluent integrates it in its original formulation [Eq. (9)] with a 5th order Runge–Kutta scheme.⁹⁴ While the particle velocity appearing in the drag contribution is treated in an implicit way, all the other terms are dealt with explicitly (which might be regarded as another difference with respect to the implementation described in Sec. III A, where all terms present in both the Navier–Stokes and the Maxey–Riley equations were handled in the framework of a time-explicit scheme).

For what concerns the evaluation of the forces acting on particles, we have selected the following options available in Fluent: the “barycentric method” to reconstruct the fluid velocity at the instantaneous positions occupied by particles (needed to evaluate the drag), and bi-linear interpolation and first-order Euler schemes to determine the material derivative DV/Dt (needed to account for the force exerted by the undisturbed flow on the particle and the virtual-added mass force).

The accuracy control implemented naturally by Fluent allows the use of different time steps for the particle motion equation and the Navier–Stokes equations. By setting the accuracy control to 10⁻⁵, we have enforced an integration time step lower than the relaxation time of the particles. With such settings, in particular, particles have been forced to perform “at least” three steps in a cell for every time step relating to the integration of the Navier–Stokes equations.

Results produced with Fluent are reported in Fig. 2 where they are directly compared with the equivalent ones yielded by the computational platform described in Secs. II and III A. The reader is directly referred to the figure caption for detailed information about the representative cases for which simulations conducted with the two platforms have been compared.

All cases have been initialized with a uniform distribution of particles. The minor visible differences between the present and the Fluent-based results should be ascribed to the slightly different initial positions considered for particles and/or to the fact that the snapshots do not relate exactly to the same times (due to the different time integration steps required by the two platforms) and/or to the different interpolation schemes used by Fluent and the present code. The very good agreement made evident by these figures (despite the significant differences in the underlying numerical strategies) should be regarded as a good indicator of the lack of implementation errors and the physical relevance of the findings. As a concluding remark for this section, we wish to highlight that other comparisons demonstrating consistency were also considered for the 3D cases (not shown).

As already explained to a certain extent in the introduction, the present work stems from the argument that a dense distribution of solid particles in a fluid undergoing thermovibrational convection may disturb the underlying transport mechanisms, thereby setting an upper limit (in terms of parameter φ for fixed particle size and density) to the effective formation of ordered particle structures.

In this section, the potential of the proposed new framework is demonstrated by considering a series of representative cases of practical complexity and relevance (i.e., suitable for the definition of real experiments to be conducted in space) and the ensuing results are presented following a specific order, which reflects the need to implement a precise analysis hierarchy. In particular, two important intertwined issues are addressed, i.e., the concurrent possibilities to influence the particle-forming mechanisms by increasing the density of particles or their overall amount, i.e., the volume fraction (both leading to an increase in the mass load).

As we will illustrate in the remainder of this section, vibrated particle-fluid systems can be extremely sensitive to these parameters and produce spontaneous “asymmetries.” Resolving these questions, therefore, also requires estimating the role played by the dimensionality of the problem (that is the number of space dimensions involved). This is the reason why some of the explorative simulations are also conducted in a two-dimensional (2D) framework (we will come back to this important concept later).

Pattern formation in these processes is described here with respect to particle motion, related structures, and the convective multicellular patterns arising as a consequence of vibrational buoyancy forces. The following values or ranges of the characteristic parameters (defined in Sec. II) are considered: Prandtl number (Pr = 6.11 corresponding to pure water at ambient temperature), non-dimensional vibration angular frequency 10³ ≤ ϖ ≤ 10⁴, Rayleigh number 10⁴ ≤ Ra_ω ≤ 10⁵, non-dimensional acceleration amplitude 10⁸ ≤ γ ≤ 5 × 10⁸, density ratio 0.1 ≤ ξ ≤ 2 and St = 5 × 10⁻⁶ (all these non-dimensional parameters reflect the space hardware and materials described in Ref. 66, namely, cubic fluid containers with size L= 1 cm, ΔT in the range 0.5–15 °C, vibrations with frequency spanning the interval 1–2.8 Hz, acceleration amplitude between 2.8 and 9.6 m/s² and particles made of glass, with various densities and diameter ≅90 μm). Though, for consistency with the companion study conducted by Lappa and Burel,⁶⁴ the Stokes number is fixed to St = 5 × 10⁻⁶, the volume fraction of particles is allowed to span a relatively wide interval, which combined with the considered set of values of ξ, leads to a wide range of mass loads (χ). We wish also to recall that, as illustrated by Refs. 56 and 57, if the density of the particles is set equal to the density of the liquid, no particle attractors exist. Particles would remain uniformly distributed in the fluid, which explains why the case ξ = 1 is not considered.

Continuing in a similar vein to Ref. 64, most conveniently, we start from a short description of the so-called time-averaged flow, i.e., the convective effects that are produced in addition to the instantaneous ones when the Gershuni number takes relatively high values and vibrations are perpendicular to the imposed temperature gradient.

In ideal conditions, that is, with no particles being dispersed in the fluid that can disturb it (see Fig. 3), this type of convection typically gives rise in rectangular cavities to the so-called quadrupolar field, i.e., a set of four rolls satisfying conditions of symmetry with respect to the horizontal and vertical axes. Such symmetry can be broken if the direction of vibrations is not perpendicular to the imposed temperature gradient and/or if the Gershuni number exceeds a given threshold. In such conditions, the original perfectly centro-symmetric distribution of rolls is taken over by a single roll pervasive throughout the physical domain with two small counter-rotating cells located in opposing corners of the cavity.^64,95,96

Apart from displaying a quadrupolar pattern, the time-averaged flow for relatively small values of Gs is also steady and relatively weak. This means that if, at the same time, the Rayleigh number [Eq. (2)] is high, the oscillatory components of the velocity field are dominant [as implicit in Eq. (3), for high Ra_ω, conditions of small Gs can be attained if the angular frequency of the vibrations is sufficiently large].

In particular, as illustrated by Ref. 56, the dominant instantaneous flow (for a square enclosure, small Gs and high Ra_ω) simply consists of a single roll occupying the entire cavity, which continuously changes its sense of rotation (from the clockwise to the counterclockwise orientation) in a consistent way with respect to the direction of the time-varying acceleration.

Without particles, indeed, the flow response to the applied vibrations is, in general, perfectly synchronous, that is, if the forcing has non-dimensional frequency ϖ, the signal measured by a generic probe put in the fluid in terms of instantaneous velocity or temperature will display exactly the same frequency (and/or frequencies that are multiples of ϖ, i.e., higher order harmonics). The signal will obviously also include a time-averaged component, which, as explained before, can take negligible or appreciable values depending on the specific regime considered. As shown, e.g., by Ref. 42, the detectable steady component in the signal is close to zero if ϖ is relatively small (in such conditions the flow is dominated by the instantaneous component of velocity). As stated before, however, the relative importance of the time-averaged flow can be weakened also for high ϖ provided Gs is sufficiently small (i.e., Ra_ω is large). As anticipated in Sec. II A, while the amplitude of the instantaneous velocity grows with Ra_ω, the magnitude of its time-averaged part increases with Gs.⁵⁶ 

Since the interplay between the non-symmetric pattern shown in Fig. 3(b) and the emerging particle structures has been already examined in a companion work,⁶⁴ in the present study we will limit ourselves to considering values of Gs ≤ O(10³) for which the time-averaged flow is weak and it always corresponds to the classical quadrupolar field depicted in Fig. 3(a).

The preliminary descriptions elaborated in Sec. IV A (dealing only with purely fluid-dynamic aspects) logically pave the way to this section, where particles are added to the fluid. As the main criterion to judge on the ability of two-way coupled particles to influence the overall system dynamics, in particular, we consider the perturbations that they can create in the flow spatiotemporal behavior. As a testbed for such analysis (without loss of generality), we consider again the classical two-dimensional square cavity (Figs. 4–10).

As an example, Fig. 4 shows the outcomes of the numerical simulations carried out neglecting the back influence of particles on fluid flow [see, e.g., Fig. 4(a) for ξ = 0.5] together with equivalent computations conducted in the frame of the two-way coupling model [Fig. 4(b)]. As outlined above, this extension to one-way coupled results is intentionally used to verify by cross-comparison the efficacy of particles in disturbing the carrier thermovibrational flow. In such a context, indeed, it is worth recalling that with the one-way coupled approach inertial particles are “merely” transported, i.e., their number has no impact on the characteristics of the carrier flow [a high number of particles being used only to reveal the perfection of the underlying attractors which behave as “templates” for their accumulation, Fig. 4(a)]. Along the same lines, we also allow [see Fig. 4(c)] the volume fraction φ to exceed the limit for which a four-way coupling approach would be required. This “license” [which we take only for the 2D results, as all the 3D results presented in the remainder of this study obviously obey the constraint φ < O(10⁻²)] is intentionally used to determine in a simplified framework the conditions in which, even by filtering out the disturbing particle-to-particle effects, the back influence of particles on fluid flow would be so large to eventually prevent the manifestation of the sought attractors.

Given these premises, some initial insights follow naturally from a comparison of Figs. 4(a) and 4(b). While Fig. 4(a) shows the classical two accumulation loops already extensively reported in previous studies (the flow being insensitive to the amount of particles^56–59), in Fig. 4(b), it can be seen that for a particle volume fraction φ = 3.6 × 10⁻³, the perfect symmetry and high-resolution appearance of the structures is taken over by a less ordered distribution. Though the existence of attracting regions is still recognizable, a series of “disturbances” can be detected. These disturbances manifest themselves in the form of “eruptions” of swarms of particles occasionally leaving the (left and right) accumulation regions [visible in Fig. 4(b) as horizontal jets, i.e., fingers of particles elongated in the direction of the applied vibrations].

Notably, by changing φ there are an uncountable number of such patterns, which differ in the degree of disorder. This is the reason why, to elucidate further the impact of this parameter, we have increased it by one order of magnitude [the limiting case with φ = 2.8 × 10⁻² shown in Fig. 4(c)]. This figure is useful as it reveals that, even if particle interaction and related effects are neglected, a too high mass load can definitely prevent the dispersed solid mass from being accumulated. This is equivalent to stating that an upper boundary in terms of φ (or χ) must be put on the region of existence of attractors potentially identifiable in the frame of one-way coupling models. In other words, this means that the “attractee” can effectively manifest in the form of well-recognizable structures of particles only if χ does not exceed a certain threshold value.

Following up on the previous point, as a fruitful alternative for obtaining additional insights into these behaviors, we have examined the velocity signals measured by a numerical probe located in a representative point (the geometrical center of the cavity). These are shown in Fig. 5.

This figure is instrumental in demonstrating that the back influence of particles on the fluid can make the resulting flow essentially non-periodic. Indeed, a series of long period disturbances can be recognized in Fig. 5(a) while Fig. 5(b) displays a much more complex frequency spectrum.

These plots can also be used for another interesting interpretation. Figure 5(c) [containing the frequency spectra related to the signals shown in Figs. 5(a) and 5(b)], indeed, demonstrates that the average amplitude of the velocity signal can significantly be increased for a fixed value of ξ on making φ larger. We can readily see the physical significance of this finding by reflecting on the transfer of momentum between solid particles and fluid. Under the effect of vibrations, particles can gain momentum directly as a result of the induced acceleration and their different density with respect to the surrounding liquid. Part of this momentum is transferred (via viscous effects, i.e., particle drag) from particles to fluid, thereby energizing the fluid flow.

Put simply, the presence of particles plays the role of an extra source of momentum for the fluid [formally represented by S_m in Eq. (15)] in addition to the classical buoyancy term (⁠ $Pr R a ω T sin ( ϖ t ) n ̂ ¯$⁠). This feedback mechanism can, therefore, be seen as the root cause for the increase in the amplitude of the signal visible in Fig. 5(b) [corresponding to the red line in Fig. 5(c)], the manifestation of which (in terms of patterning behavior) is the recognizable strong asymmetry in the distribution of particles [Fig. 4(c)].

Another meaningful way to think about the back influence of particles is to consider the effect they can exert on the fluid from a time-averaged standpoint. This aspect is quantitatively and qualitatively substantiated by Fig. 6 for φ = 3.6 × 10⁻³.

Though the classical four cells typical of the quadrupolar field can somehow still be identified, their morphology is highly distorted. Moreover, topological defects, such as cracks, cusps, and discontinuities, start to affect the streamlines. They witness directly the ability of particles to exchange momentum with the fluid; in Fig. 6, it can also be seen that the quadrupolar field for the considered value of Gs is no longer time-independent, i.e., it changes with time.

Another related key observation concerns the strength of the time-averaged convection, which tends to be higher when particles are present (ψ_max = 8 × 10⁻² in place of ψ_max = 1 × 10⁻² without particles). This provides additional independent confirmation that on average there is a net transfer of momentum from particles to the fluid.

The situation dramatically changes when φ is increased by one order of magnitude (φ = 2.8 × 10⁻² in Fig. 7). The quadrupolar field is no longer recognizable. The time-averaged convection is characterized by the continuous nucleation of cells that with time coalesce or split resulting in a more or less chaotic behavior (Fig. 7).

Continuing with the description of the numerical results, Fig. 8 may be regarded as the analog of Fig. 4 for a value of the density ratio larger than one (ξ = 1.5), i.e., for particles denser than the fluid. Direct comparison with the situation where ξ = 0.5 indicates that the disturbing effect exerted by particles for a fixed value of the particle volume fraction φ can be enhanced by increasing the density ratio. This obvious finding simply reflects the corresponding increase in the mass load χ = φξ. Along these lines, what is emphasized in Fig. 9 is still the inter-relatedness of the disorder detectable in the particle distribution and the velocity signals (and related frequency spectrum).

By taking a look at the vertical axes, the reader will easily realize that the fluid motion is yet energized by the interplay with the dispersed solid mass. The frequency spectrum, however, is more involved, which indicates that the disturbances preventing the flow from taking a regular behavior are more complex and intense [compare the red lines in Figs. 9(c) and 5(c)]. This trend is further confirmed by the patterns reported in Figs. 10 and 11 for φ = 3.6 × 10⁻³ and φ = 2.8 × 10⁻², respectively.

The simplification to two dimensions at the basis of all the results presented in Secs. IV A and IV B may seem natural at first sight because of the intrinsic symmetries of the system sketched in Fig. 1. Nevertheless, it should expressly be pointed out that, since the two-way coupled evolution of a non-isothermal fluid and dispersed solid particles under the effect of vibrations has not been investigated until now,^56–59 there is no reason to think a priori that this interplay should not hide completely unknown mechanisms somehow involving the spanwise direction. This is the reason why we have intentionally based the present work on an analysis hierarchy implemented through a characteristic path of progression from 2D to 3D simulations. First two-dimensional simulations have been executed (where any processes that depend on the details of the third direction are excluded, Sec. IV B); then, three-dimensional simulations have been conducted (this section) to reveal heretofore unseen 3D effects via intentional cross-comparison with the equivalent 2D ones.

Following the same approach undertaken in Sec. IV B, in particular, in this section, first we discuss cases for which particles are lighter than the external fluid and then we focus on situations with ξ > 1. Along these lines, the outcomes of the 3D simulations for the same circumstances already considered in Fig. 4(b) are reported in Fig. 12 (3D snapshot of particle distribution) and Fig. 13 (related time evolution).

As a fleeting glimpse into Fig. 12 would confirm, symmetry breaking phenomena with respect to the spanwise direction can effectively occur when particles have an additional spatial degree of freedom (i.e., they can move freely along x). This behavior becomes very evident if the distribution of particles is inspected by taking a look at it along an axis perpendicular to the xz [Fig. 12(b)] and xy [Fig. 12(c)] planes, respectively. Particles do not remain uniformly distributed along x (as the assumption of 2D behavior would require); rather, they display some preferred planes of accumulation. More specifically, three disjoint main planes supporting clustering dynamics can be distinguished in Figs. 12(b) and 12(c), one approximately located at the center of the cubic enclosure and two located in proximity to the walls delimiting it along the x direction.

The 3D self-organization of the particles also implies remarkable changes in the pattern that can be seen considering the projection of the solid mass distribution in the yz plane [compare e.g., Fig. 12(a) with Fig. 4(b)]. Particles do still tend to cluster forming structures that an external observer (looking at the fluid domain along the x axis) would recognize. However, the morphology of these aggregates is quite different; the areas where solid matter is concentrated have an irregular triangular shape [in place of the elongated elliptical accumulation regions visible in Fig. 4(b)].

To put these results in a broader perspective, in the following we will attempt an interpretation of such variations in the morphology of the particle structures on the basis of the “bifurcation theory.” Indeed, it is widely recognized in fluid dynamics (and more, in general, in the field of non-linear systems) that the loss of symmetry implies the existence of a new solution that bifurcates from a preexisting one due to the selection and ensuing amplification of disturbances (in this regard the solution shown in Fig. 4(b) could be regarded as the “basic 2D state” from which the 3D state shown in Fig. 12 originates as time increases).

Although the analogy with the classical concept of bifurcation in fluid dynamics should effectively be considered as a possible key for the interpretation of these results, the identification of the related physical “disturbances,” that is, the cause-and-effect mechanisms responsible for the symmetry breaking process, however, is not as straightforward as one would assume. Since for the considered value of the Gershuni number no 3D instability occurs in pure thermovibrational flow (i.e., when no particles are considered) or when the back influence of the dispersed solid matter on the carrier flow is not taken into account (one-way coupling), we argue that the sought mechanism implicitly relies on the two-way coupling between fluid flow and particle motion, that is, the 3D disturbances take their energy from such an interplay.

Interestingly, this may be regarded as another typical example of systems where separated entities can produce collective behaviors more complex than those of the individual components (see, e.g., Ref. 97). As revealed by the present 3D simulations, two-way coupled particles can provide these vibrated dispersions with the required internal feedback loops by which they can develop their own capacity for transformation, requiring only the right conditions for activation. As already described to a certain extent in Sec. IV B, these processes (able to feed information back into the system where it is iterated, or used multiplicatively) are established in terms of momentum exchange.

The continuous exchange of momentum between fluid and particles enabled by particle drag and inertial effects can make these processes extremely sensitive to their (even though very small) internal variations allowing the amplification of initially extremely small disturbances. A clear separation should therefore be considered between these results and those previously obtained in the frame of the one-way approach by Lappa and Burel,⁶⁴ where the observed 3D phenomena were essentially a consequence of the three-dimensional nature of the time-averaged field established inside the cavity for relatively high values of the Gershuni number [Gs = O(10⁴)].

A snapshot of the numerically computed disturbances for the present case (that we represent in terms of the time-averaged fluid velocity component along the spanwise direction) can be seen in Fig. 12(d). The highly localized nature of these fluid-dynamic perturbations further supports our interpretation about the role played by the exchange of momentum between fluid and particles in inducing the 3D instability (indeed, the regions of the physical space where the fluid velocity spanwise component is significant closely mimic the distribution of particles).

Further pursuing the possible analogy with the typical mechanisms of amplification of disturbances envisaged by the linear stability analysis (where the growth of perturbations is assumed to follow an exponential law, see, e.g., Ref. 96), Q⁺ has been reported in Fig. 13(a) using a logarithmic scale. As the reader will easily realize, using such a scale the perturbations indeed grow following an apparently linear law, which confirms the exponential behavior of the disturbance growth process. This figure is useful as it also shows that after an initial transient time, disturbances saturate their amplitude and a new stable state is attained, which is consistent with the theory of bifurcations.

To do justice to these interesting findings, Sec. IV D is entirely devoted to a characterization of this bifurcation and identification of the conditions that make it possible.

Starting from the arguments elaborated in Sec. IV C about the essentially inertial nature of the instability producing transition to 3D flow, the next natural step toward a complete characterization of this phenomenon would require proper assessment of the sensitivity of the symmetry breaking process to typical problem parameters, especially those that can exert an immediate impact on the momentum exchange between the two phases [namely, the particle volume fraction (φ), the acceleration amplitude (γ), the density ratio (ξ), and the Rayleigh number (Ra_ω)].

As shown in Fig. 14(a), a decrease in the volume fraction of dispersed solid mass for a fixed value of ξ can change the time required for disturbance saturation, which provides additional important clues regarding the essentially inertial nature of the instability.

Most interestingly, making φ smaller can even lead to conditions for which the emerging structure recovers the perfect geometrical morphology of the underlying attractors determined in the framework of one-way coupled numerical studies [Fig. 14(b) to be compared with Fig. 8(a)]. This figure also quantitatively substantiates the scarce role played in this case by the time-averaged fluid velocity component in the spanwise direction (it being limited in strength and present only in very localized regions).

Notably, similar effects can be obtained by properly tuning the acceleration amplitude. In this regard, Fig. 15 is extremely instructive as it proves that other situations exist for which the 3D disturbances are not excited and the pattern can retain its essentially two-dimensional configuration. These figures refer to the same value of the density ratio and particle volume fraction already considered for Fig. 12 (ξ = 0.5 and φ = 3.6 × 10⁻³), i.e., the same mass load (χ = φξ), but smaller values of γ and ϖ (γ = 1 × 10⁸, ϖ = 5 × 10³).

We did not detect any symmetry breaking process for such conditions (the reader being directly referred to Figs. 15(a) and 15(b) for evidence of the essentially 2D behavior, witnessed by the order of magnitude of Q⁺ and of the fluid velocity in the center of the cavity, respectively). Although a 3D disturbance exists [Fig. 15(a)], its amplitude and magnitude are relatively limited [O(10⁻²)]; moreover, it displays no tendency to be amplified as time increases even if the simulation is prolonged up to 10 times the thermally diffusive reference time (for a cubic cavity having a size of 1 cm, this would correspond to a physical time of approximately two hours).

Along these lines, further understanding of the specific equilibrium conditions attained in this case can be gained by considering again simulations conducted under the constraint of two-dimensional flow (Figs. 16–18).

These results are meaningful as they can be used yet to easily get quantitative information on the inertial mechanisms at work in the considered process. As a quick look at Fig. 17(a) would immediately confirm, the response of the fluid–particle mixture to the applied forcing is very regular for φ = 3.6 × 10⁻³. For the limiting condition represented by φ = 2.8 × 10⁻², though some spikes can occasionally be seen in the signal, the behavior is still essentially time-periodic [Fig. 17(b)]. The conceptual ingredient needed to understand these results lies in considering the impact of the parameter γ on particle acceleration [via Eq. (17)] and on the associated momentum exchange term at the right hand side of the Navier–Stokes equations [represented by Eq. (20)]. By causing a decrease in the momentum transferred between the two phases, a decrease in γ is beneficial in terms of stability behavior because it limits the disturbances that particles can induce in the underlying fluid flow. Similar concepts also apply to the time-averaged field. The quadrupolar field visible in Fig. 18(a) can clearly be recognized even if the particle volume fraction is largely increased [Fig. 18(b)]. These observations align with the trend visible in terms of patterning behavior in Fig. 16.

By limiting the energy that the system can use for the amplification of disturbances, a decrease in γ can make the overall system stable, this argument being fully supported by the 3D simulations [Fig. 15(c)].

There is also another interesting extension that attaches to this interpretation. An increase in ξ should, indeed, exert an opposite influence on the dynamics, i.e., it should promote the onset and growth of disturbances for a fixed value of γ (that is why in the remainder of this section we also probe the role of this factor). As an example, the outcomes of the 2D simulations for Ra = 10⁵, γ = 1 × 10⁸, ϖ = 5 × 10³ and ξ = 2 are summarized in Figs. 19 and 20.

As evident in these figures, for the same value of the acceleration amplitude, vibration frequency, Rayleigh number, and particle volume fraction (φ = 3.6 × 10⁻³), but a higher value of ξ [corresponding to four times the mass load relating to Fig. 16(b)], the feedback loop by which disturbances can be excited starts to be effective. This is consistent with the presence of jets of particles being erupted from the structures in Fig. 19(b) and by the alterations undergone by the classical quadrupolar field (see Fig. 20, interestingly, in these circumstances, the back influence of particles on fluid flow causes quite regular periodic coalescence and separation of time-averaged rolls located in opposing corners).

The increased interference that dispersed mass has with the spatiotemporal behavior of the fluid-dynamic field is confirmed by the 3D simulations (Figs. 21 and 22). As shown by Fig. 21, an initial transient phase exists where disturbances exhibit a more or less two-dimensional nature, resulting in the eruption of jets of particles in a direction parallel to the imposed vibrations [Fig. 21(a), left panel] already seen in the 2D simulations. At a certain stage, however, the “curse” of dimensionality starts to affect the dynamics and 3D disturbances tend to be amplified. The related stages of evolution (yet collected in Fig. 21) are extremely interesting.

In particular, as made evident by the views along a direction perpendicular to the xz plane (central panels of Fig. 21), the set/distribution of particles, initially uniformly spread along the spanwise direction, with time tends to split into two main distinct “swarms.” At first, the two swarms are very close each other [giving the illusion of a single spatially extended bunch of floating particles occupying the center of the cavity, Fig. 21(b), central panel]. As time passes, however, the two groups of particles come apart (their separation along the x direction tends to increases). This process is accompanied by a reorganization of the dispersed mass, which progressively collapses into two thin accumulation regions stretched in the direction of vibrations, i.e., the y axis [hereafter simply referred to as “accumulation planes,” see Fig. 21(c), central panel]. The distance between these two planes keeps increasing with time. As a result, the two thin accumulation regions are progressively displaced toward the two walls delimiting the cavity along the x direction. Finally (when the amount of clear fluid located between the thin accumulation regions and these walls becomes relatively small, i.e., when the accumulation planes occupy a position relatively close to the end wall), the particle-dense regions lose their elongated shape and are turned into two “blobs” located in opposite corners of the cubic cavity [as seen in the xy plane, Fig. 21(d), right panel for t = 10]. Moreover, an observer looking at the fluid container along the x direction at this time would see triangular-like structures resembling those already found for γ = 5 × 10⁸, ϖ = 10⁴, ξ = 0.5, φ = 3.6 × 10⁻³ [Fig. 21(d), left panel].

Remarkably, this entire process can directly be put in relation with the time-averaged component of fluid velocity in the spanwise direction. As witnessed by Fig. 21(e), indeed, it is the time-averaged component induced by the transfer of momentum from particles to fluid that exerts a compressive action on the particles forcing them to cluster into regions of limited extension (small thickness) along the spanwise direction and determining the progressive displacement of these particle-dense loci toward the external walls.

Having completed a sketch of the potential impact of the acceleration amplitude (γ), the particle volume fraction (φ), the particle-to-fluid density ratio (ξ) (and therefore the mass load χ = ξφ), we now turn to assessing the role of the last influential factor, that is, the Rayleigh number (Ra_ω).

While the previous parameters were all directly affecting solid matter dispersed in the fluid, this number will be exerting an influence on the dynamics following a different cause-and-effect connection, i.e., due to its direct relationship with the motion of fluid. A decrease in this parameter will clearly produce a mitigation of the fluid flow strength and therefore weakening of the inertial forces at play in the considered dynamics.

For the same values of ξ, φ, and γ used for the simulations reported in Fig. 21 (namely, ξ = 2, φ = 3.6 × 10⁻³ and γ = 1 × 10⁸), consideration of smaller values of Rayleigh number and angular frequency (Ra_ω = 10⁴ and ϖ = 10³, respectively) can have (as expected) remarkable beneficial effects on the regularity and perfection of the emerging structures (Figs. 23 and 24). Although some 3D disturbances can still grow (Fig. 23), particles can successfully reveal the perfection of the underlying attractors (Fig. 24).

Interestingly, the particle distribution reorganization process follows a slightly different evolution path in this situation. Like the case discussed before (γ = 1 × 10⁸, Ra_ω = 10⁵, ϖ = 5 × 10³, ξ = 2.0 and φ = 3.6 × 10⁻³), as time passes, the initially continuous distribution of particles along the spanwise direction is broken into two different patches [the entire process may give the observer the illusion of a “cellular mitosis,” as evident in Fig. 24(b), central panel]. The formation of disjoint accumulation regions spanning the entire extension of the cavity in the direction of vibrations, however, does not occur in this case. The particle-dense regions are not stretched in the y direction. Rather, as a result of the entire process, two opposing blobs of particles emerge directly. As shown in Fig. 24(c) (right panel), each blob tends to occupy the corner between one of the walls limiting the cavity in the spanwise direction and another wall limiting it in the direction of vibrations. These two compact accumulation regions are mirror symmetric with respect to the diagonal direction x = y.

As a concluding remark for this section, Figs. 25 and 26 witness that, in line with the various interpretations elaborated before, disturbances can be further mitigated by decreasing the density ratio. For γ = 10⁸, ϖ = 10³, Ra_ω = 10⁴, and ξ = 0.5, particle structures can again manifest as geometrically perfect entities mirroring in the physical space the perfection of the corresponding attractors existing in the one-way coupled space of phases.

In earlier studies, the relationship between “particle attractee” and thermovibrational flow has been investigated under the assumption of one-way coupling (which means that attractors have been determined as abstract entities playing the role of “templates” for the accumulation of passively transported inertial particles). Stripped to its basics, the present research scheme has envisioned a parametric investigation aimed to identify the conditions where these ideal entities can still manifest in the form of recognizable particle structures (with a more or less high degree of perfection) if the assumption of dilute dispersion ceases to be valid.

At present, there are no published papers examining this issue; hence, this is the first time that detailed data on the related phenomenology (morphology and range of existence of the structures, 3D velocity disturbance growth rates and eventually complete transition to chaotic states) have been presented.

Starting from the realization that earlier models used to address this problem had to be modified to allow for realistic particle structures to exist, the first objective of this work was to pursue modeling of a more complete and sophisticated theoretical and numerical framework. This has been achieved by properly coupling Eulerian (for fluid flow) and Lagrangian (for particle tracking) approaches in such a way to allow realistic inter-phase feedback effects.

In addition to purely mathematical and numerical aspects, we have implemented a specific analysis hierarchy by which the problem has partially been tackled under the unphysical constraint of two-dimensional flow. These simulations have successfully been used as a workhorse to test the efficacy of two-way coupled dispersed particles in disturbing the carrier fluid flow. In particular, under this constraint, circumstances have been revealed where the back influence of dispersed solid mass on the carrier flow can even overshadow the underlying attractors (predicted on the basis of one-way coupling studies) and lead to completely chaotic particle dynamics.

Remarkably, the 2D study has also been instrumental in unraveling interwoven or overshadowed processes, which would have otherwise been difficult to interpret, such as the unexpected rupture of symmetry occurring in the spanwise direction.

Using a synergetic combination of 2D and 3D computations, we have conducted an analysis of the triadic relationship among the intensity of fluid convection, the amplitude of acceleration and the frequency of the vibrations. Notably, this interplay takes place within an interlocking ensemble of mutual interferences, where the primary relationship between the fluid flow and its ability to transport the particles (depending on their geometrical and physical properties, i.e., the Stokes number and density ratio) is mediated by the secondary influence that particles have on the fluid flow itself (where also the particle volume fraction starts to play a significant role), as well as the tertiary influence that the flow disturbed in such a way, in turn, has on them.

In turn, the tertiary influence encompasses a new category of three-dimensional bifurcations. These are responsible for the emergence of a velocity component along the spanwise direction. This kind of instability is essentially driven by the momentum exchange between fluid and particles and the ability of the latter to energize the carrier flow. In this regard, this process may be regarded as a new (inertial) mechanism able to produce symmetry breaking phenomena, which have nothing to do with the equivalent ones induced by the time-averaged flow when the Gershuni number is relatively high (effects which would be present even in the complete absence of dispersed solid matter, see, e.g., Ref. 64).

In order to obtain an indirect proof of such conceptual scenario, we have intentionally kept the value of the Gershuni number relatively small and changed the parameters potentially influencing the inter-phase exchange of momentum, namely, the Rayleigh number, the acceleration amplitude and the particle density and volume fraction. Such analysis has led to the conclusion that while the 3D instability can be promoted by an increase in the density ratio and/or the particle volume fraction, vice versa a decrease in the acceleration amplitude and/or Rayleigh number can stabilize the carrier flow.

Discerning the role played by the effective percentage of mass dispersed in the fluid has contributed to our understanding of these dynamics and to the capacity to optimize future experiments planned for execution onboard the International Space Station. Perhaps most important of all has been the identification of new routes to as-yet unknown phenomena. We have shown that, as a result of the above-mentioned 3D instabilities, new types of structures can be formed, which differ from those reported in earlier work on the subject due to their prevailing triangular morphology (as seen by an observer taking a look at the flow along a direction perpendicular to the plane of vibrations and imposed temperature difference) and owing to their ability to collapse in planes of accumulation or “blobs” well separated along the spanwise direction.

Future studies shall be devoted to the investigation of all these dynamics when two different mechanisms are at play at the same time and determine 3D flow, namely, the phenomena of purely fluid-dynamic nature that are produced when the Gershuni number takes relatively high values (able to break the symmetry of the quadrupolar field and produce significant values of the time-averaged velocity component in the spanwise direction⁶⁴) and the 3D bifurcations driven by fluid–particle (two-way) momentum transfer revealed by the present work.

This work has been supported by the UK 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 author would like to thank Dr. Thomas Burel for the additional simulations presented in Sec. III B conducted using the ANSYS Fluent computational platform and Alessio Boaro for computing the spectra shown in Figs. 5(c) and 9(c).

The authors have no conflicts to disclose.

The data that support the findings of this study are openly available in Pure at https://doi.org/10.15129/b70c149c-04f2-47a6-b23d-50baa35da491, Ref. 99.