Vibrationally driven particle formations in fluid systems with bimodal thermal inhomogeneities

Particle “attractors” in liquid flow, i.e., loci where solid particles tend to accumulate as time passes, are essentially a consequence of a mismatch in the fundamental nature of the liquid motion and that of the transported particles. More precisely, while the former is essentially incompressible, i.e., from a purely mathematical standpoint, it satisfies the identity “divergence of velocity equal to zero” (corresponding from a physical point of view to the conservation of volumes), the latter can violate this constraint thereby allowing the space delimited by different particles to shrink or grow. Notably, in the first case, this represents the necessary (but not sufficient) pre-requisite allowing the existence of particle attractors in the fluid flow. A swarm of particles can behave as a kind of “compressible medium” essentially because, owing to inertial and drag effects induced by the finite size and mass of the particles, their trajectories are not forced to mimic exactly those of the carrier fluid flow. Small departures are produced and accumulate with time, thereby allowing a swarm of particles to be compressed or expanded in time.

This category has been placed in a precise theoretical context in a series of works (see, e.g., Lyubimov et al.,^1,2 Haller and Sapsis,^3,4 Sapsis and Haller,⁵ and Balasuriya et al.⁶) where precise theorems have been formulated to characterize these dynamics. Most interestingly, it has been shown that additional ingredients required for the effective manifestation in the physical space of these attractors (which exist in the space of phases as mathematical items, but show up in the real space as particle aggregates) is the existence of well-defined “repetitive” behaviors in space and/or time in the considered carrier flow. Such a remarkable theoretical foundation has led to various lines of inquiry running in parallel in the literature and differing with regard to the nature of the repetitive behaviors mentioned above, which allow particle clustering. These are typically represented by closed “tubes” in space, which serve as templates for the accumulation of particles or cyclic phenomena in time, by which the small departure of the particle trajectories from the fluid streamlines can accumulate in time and grow, thereby allowing particles relative distances to change.

Here, we limit ourselves to mentioning the particle attractors, which have been discovered in liquid bridges featuring time-periodic (unstable) Marangoni flow^7–16 and those revealed by other studies based on numerical simulations and experiments where the cyclic nature of the fluid flow has been artificially induced by imposing vibrations on a non-isothermal flow.^17–23 Notably, these phenomena require the considered system to be finite in space as the interaction of particles with the boundary somehow plays a relevant role in supporting their aggregation (a concept that applies to both thermocapillary and thermovibrationally induced particle structures).

Along these lines, the present study may be regarded as another example pertaining to the line of research where externally used vibrations are used to drive the phenomena of interest.^17,22 This problem is somehow unique given the additional flexibility brought in by the several possible degrees of freedom related to the possible thermal boundary conditions and direction of the driving force with respect to other cases (e.g., the above-mentioned liquid-bridge problem). Although such phenomena were originally discovered for a differentially heated cavity with a uniformly heated wall and an opposing uniformly cooled wall,^17,18 and later studies could show that with such a configuration a variety of particle structure shapes can be obtained depending on the relative direction of the imposed temperature gradient and vibrations,^19,20 there is no specific reason for which a uniform temperature along the walls should be considered.

This led some investigators^24,25 to explore circumstances where thermal inhomogeneities were deliberately introduced on otherwise uniformly cooled or heated walls. In turn, this led to the realization that the considered problem has many facets, which could be exploited to support current trends toward the identification of novel manufacturing methods. Many materials used in our everyday life pass through a fluid state (before being solid), which consists of fine particles or droplets (the minority phase) distributed inside an external (liquid) matrix (the so-called majority phase). Having the ability to control the distribution of solid particles with no need to physical contact or manipulate them may lead, in principle, to a variety of new materials with heretofore unseen properties.^26,27

Lappa^17,19 reported on the ability of vibrations to produce structures with the morphology of the quadrics surfaces (ellipsoids, cylindrical surfaces, paraboloids, conical surfaces, etc.) in uniformly heated/cooled cavities and clarified that such structures are always obtained as a couple of formations with a certain symmetry coexisting in the considered containers. Lappa and Burel²⁸ investigated the dependence of the structures formation time and characteristic size on a number of influential parameters, these being the nondimensional vibration acceleration amplitude, frequency, Rayleigh number, Stokes number, and particle-to-fluid density ratio. Accordingly, a general correlation was introduced to predict the formation time in various circumstances. Lappa²⁹ assessed the potential impact of the number of particles, i.e., the more or less dilute or concentrated nature of the considered dispersion on the effective particle structures formation. They found that on increasing the number of particles beyond a certain threshold, which depends on the nondimensional vibration acceleration amplitude and other parameters, a discrepancy is produced between the theoretical attractors revealed by one-way coupled numerical simulations and the effective aggregations showing up in the physical space. The latter tend to be less regular, i.e., to display a corrugation that increases with the particle concentration and/or density. Moreover, the back influence of the particles on the fluid flow may cause symmetry breaking instabilities, which would not be possible in the absence of particles.

Most recently, Crewdson et al.²⁴ and Santhosh and Lappa²⁵ reported that the multiplicity of structures, i.e., the number of disjoint clusters showing up in the physical space can be increased if the (otherwise uniformly) heated or cooled walls are perturbed via the introduction of a central thermal spot having opposite sign (i.e., a cold spot at the center of a hot wall and vice versa). More specifically, Crewdson et al.²⁴ concentrated on two-dimensional cavities, while Santhosh and Lappa²⁵ extended the study to three-dimensional (3D) systems showing that not only non-uniform thermal boundary conditions with central spots of varying sizes can increase the number of independent formations but they can also support new accumulation mechanisms, which are not possible in uniformly heated and cooled cavities.

In the present study, we probe another influential parameter, that is, the number of thermally controlled walls, by increasing them from two to four walls as shown in Fig. 1. In order to reduce the complexity of the problem, following the same approach implemented by Refs. 24 and 25, the subject is approached in the framework of one-way coupled simulations, thereby filtering out possible flow instabilities and allowing an analysis of the theoretical attractors, which serve as templates for the accumulation of particles.

This analysis is deliberately kept short as it should be regarded as an attempt to elucidate the effect of an influential parameter only and, in this regard, it should be seen as a focused extension Ref. 25.

As already discussed to a certain extent in the introduction, the fundamental mechanisms driving the phenomena of interest are essentially the same examined by Santhosh and Lappa.²⁵ Here, however, the degree of complexity of the thermal boundary conditions is increased. Although the same cubic cavity already considered in many earlier studies pertaining to this line of inquiry is assumed, the number of thermally controlled faces is increased from 2 to 4 (see Fig. 1). Only two walls retain an adiabatic behavior, while the other four boundaries consist of two adjacent thermally heated walls, each with a central cold spot of variable size and other two (opposed) adjacent walls that have a cold temperature, with the exception of the two related central heated spots. Clearly, this situation is not equivalent to that where walls with similar temperatures are faced opposite to each other, which will be addressed in a future study (together with the companion problem where the thermal inhomogeneities are located at the corners of the cavity rather than in the wall centers). In all cases, no gravity is present, i.e., the considered phenomena are assumed to occur in microgravity conditions.

In analogy with the existing literature on the subject,^17–23 the problem is approached by tracking the fluid and the related particle approach through well-defined and already validated theoretical frameworks, namely the Navier–Stokes equations for the liquid phase and the Maxey–Riley equation for the dispersed solid matter.

The governing equations for the fluid phase have been integrated in time and in space in the framework of a projection method. With this approach, the original system of equations is manipulated mathematically in order to obtain an alternate set of partial differential equations that are better suited for the implementation of a time-marching algorithm.³³ 

For the present problem, we have selected the required grids using the criteria defined by Ref. 18, and used standard central differences to discretize accordingly both convective and diffusive terms appearing in Eq. (18). The values of the fluid velocity in the different grid points have been interpolated to determine the terms appearing at the right-hand side of Eq. (7) at the instantaneous particle position. More specifically, we have exploited simple linear interpolations schemes and integrated the particle-tracking Lagrangian equation for each particle by means of a fifth-order accurate Runge–Kutta algorithm. For each 3D simulation, 6.4 × 10⁴ particles have been tracked.

As a natural extension of the earlier investigation,²⁵ here we concentrate on the same Prandtl number (Pr = 6.1 representative of pure water at ambient temperature). In the light of the general correlation introduced by Lappa and Burel²⁸ to predict the particle structure time as a function of various influential parameters, in the present study a relatively small value of the angular frequency is considered (i.e., Ω = 10³ in order to minimize the transient period preceding the formation of well-defined structures and limit the otherwise prohibitive required computational time). For consistency with Ref. 25, moreover, we assume (γ = O(10⁸), Ra_ω = 10⁴, St = 5 × 10⁻⁶, and a particle-to-liquid density ratio ξ = 2. A distribution of evenly spaced particles is used as initial condition for all the simulations.

For the convenience of the reader, we wish to recall that under the effect of fluid drag and vibrationally induced acceleration, particles demix gradually from the fluid (where they were initially evenly dispersed) and accumulate forming the aforementioned recognizable structures; however, once such structures are formed, their shape stops changing (their morphology becomes “frozen”) and particle motion reduces to a periodic back and forth motion of the entire particle pattern along the direction of vibrations. A detailed explanation of this mechanism and of the particle dynamics in the transient stage preceding the formation of structures can be found in Ref. 17.

Yet, to support reader's understanding of the considered dynamics, here we start with the fundamental situation where the imposed vibrations are directed exactly along the y direction (the angle ϕ is equal to zero, see again Fig. 1), i.e., they are perpendicular to the only two existing adiabatic walls. These circumstances correspond to the case in which the vibrations are essentially orthogonal to the dominant temperature gradients established inside the cavity. The outcomes of the related simulations are shown in Figs. 2–7 for different values of the spot size l_s. Each of these figures consists of four panels, i.e., views showing the structures and the temperature distribution in the xy, xz, and yz midplanes and a 3D view of the structures, respectively.

In particular, the first figure of the sequence (Fig. 2) refers to the case l_s = 0.1.

Direct comparison of this figure with the equivalent one included in the earlier numerical study by Santhosh and Lappa²⁵ (Fig. 14 in that work) is instrumental in showing that the addition of two other non-uniformly heated walls has significant effects on the emerging particle formations.

In place of two almost cylindrical aggregates aligned with the z axis, complex structures are obtained. Rather than being oriented along the z axis, these align with the x = z direction and extend from one of the corners to the opposite one [Fig. 2(b)]. Furthermore, each structure gives the observer the illusion of being made up of two coaxial conical surfaces having a common base and pointing toward opposite directions (their apex being close to one of the cavity corners). This makes the overall multiplicity (N) of the pattern equal to 4, that is two times the equivalent one (N = 2) that Santhosh and Lappa²⁵ obtained in equivalent conditions with two thermally active walls only. As another characteristic parameter, we also refer to the angle φ shown in Fig. 2(b), that is, the projection of the apex solid angle in the xz plane. In the present work, it has been determined through a Matlab script relying on trigonometric formulas and particle co-ordinates pertaining to the external boundary of the structure in proximity to the corner of the cavity.

Notably, an increase in the spot size (l_s = 0.3 in Fig. 3) leads to a shrinkage in the solid angle related to the apex of each conical surface and the related projection φ. Overall, this also causes a decrease in the transverse size of such structures. However, no change occurs in the multiplicity (N = 4 as for l_s = 0.1).

Most interestingly, if l_s is increased to 0.4, important topological changes show up. The transverse size of aforementioned conical structures is reduced even more [such structures are hardly recognizable in the 3D view, see in particular, Figs. 4(a) and 4(d)]. Moreover, much more extended formations manifest in the center of the domain. Unlike the conical structures, they have a much more rounded topology and occupy vast sub-regions of the cubic cavity.

These become even more evident for l_s = 0.5 (Fig. 5), where their presence is witnessed by the approximately elliptical shape of their projections in the xy and xz planes [Figs. 5(a) and 5(c), respectively].

Interestingly, for values of the spot size larger than 0.5, the configuration with conical structures is recovered and maintained for larger l_s (see, e.g., Fig. 6 for l_s = 0.7). However, as a distinguishing mark, the related solid angle is much larger [the related projection in the xz plane is 90°, see Fig. 6(b)].

As yet evident in Fig. 6(b), in particular, the 3D surface generatrix [see Fig. 6(b)] may be regarded as the combination of a segment of limited extent parallel to the solid wall and a curved line forming with it a cusp point (indicated by the black arrow).

On increasing l_s, however, the extension of this segment becomes smaller until it disappears completely for l_s = 0.9 [Fig. 7(b)]. Accordingly, in this case a configuration formally similar to that shown in Fig. 3 is recovered.

For the convenience of the reader, such trends have been summarized in Fig. 8. The significance of this figure essentially resides in its ability to make evident the non-monotonic dependence of both the structure maximum transverse size (quantitatively measured considering its projection in the xz plane) and the aforementioned apex angle on the spot size l_s.

Moreover, notably, what stands out from Fig. 8(a) is that the considered 3D dynamics seems to display a notable analogy with those originally discussed by Santhosh and Lappa²⁵ in the 2D (square cavity) case (the interested reader being referred to Fig. 6 in their work).

The present findings, indeed, reveal a strong variation in the structure size which increases from a relatively small value to a larger one as the parameter l_s exceeds a given threshold (l_s ≅ 0.3), which may be regarded once again (as in Ref. 25) as a clue or hint for a different patterning mechanisms at work (the reader being referred to Sec. IV B for additional details about this specific aspect). Another analogy relates to the tendency of the structure characteristic size to attain similar values at the extremes of the considered range of spot sizes (i.e., for l_s ≅ 0 and l_s ≅ 1).

An additional key observation stemming from Fig. 8(a) concerns the strong decrease undergone by the structure transverse size as γ is increased, which confirms the tendency of these structures to become progressively more compact for larger acceleration amplitudes. Interestingly, a non-negligible displacement of the graph maximum toward smaller values of l_s can also be discerned, which indicates that for a larger γ, a smaller value of the spot size is needed to maximize the structure transverse extension.

Other interesting changes produced by a modification of γ or l_s are quantitatively substantiated in Fig. 8(b) in terms of the apex subtended angle. While for γ = 5 × 10⁷, this angle is always close to 90°, which may be regarded as a consequence of the much larger spatial extension of the structure and the ensuing tendency of its boundary to align with the solid walls in proximity to the cavity corners, for γ = 1 × 10⁸, a completely different trend is obtained. For this value of the acceleration amplitude, with the exception of a central interval, where as shown in Fig. 8(a) the size of the structure is comparable to that obtained for γ = 5 × 10⁷, external ranges exist where this angle attains relatively small values. This angle can be as small as 18° for l_s = 0.3 (corresponding to the particle distribution shown in Fig. 3 and the related “cusp” visible at the extremities of the conical structures). This trend is akin to the behavior already discussed for the structure size [Fig. 8(a)]. Put simply, if starting from this situation (l_s ≅ 0.3) l_s is increased, a sudden jump in the value of the angle takes place and this occurs in conjunction with the more rounded topology of the relatively large formations visible in Figs. 4 and 5 (for which an interpretation is provided in Sec. IV B).

As a concluding remark for this section, we wish to highlight that another jump (this time from an angle ≅90° to a much smaller value, i.e., ≅40°) manifests again as l_s becomes larger than 0.8, which may yet be regarded as a significant variation in the cause-and-effect relationships driving the particle self-organization process (more precisely as a return to the same formation mechanism at play for l_s ≤0.3 as illustrated in detail in Sec. IV B).

Building on the arguments originally provided by Lappa¹⁷ and much more recently by Ref. 25, in this section an interpretation for the observed trends and the related changes in the patterning behaviors is sought through direct comparison of the emerging structure morphology and the “topology” of the underlying flow field (multicellular nature and related evolution in time).

In this regard, it should be briefly recalled that two main mechanisms driving particle accumulation exist in the considered conditions. One, originally identified in the limiting condition in which no thermal spots are present (uniformly heated and cooled walls¹⁷) requires the particles to periodically hit the solid walls (owing to their different density with respect to the liquid and the force that vibrations produce accordingly) and being accumulated in straight particle-dense lines parallel to the wall itself. These are then progressively bent by the flow, thereby producing compact rounded formations. Another mechanism, reported more recently,²⁵ does not require the particles to hit any boundary. However, multicellular flow structures are needed, by which some particles can be trapped in a central region. These are typically produced when thermal inhomogeneities are considered, i.e., thermal spots are present along the sides of the system. The mark distinguishing the two types of particle structures enabled accordingly essentially relates to their external boundary. While in the former case the formations display a well-defined (particle-dense) boundary separating them from the clear fluid outside, in the latter, this boundary does not exist.²⁵ 

Notably, Fig. 9 for the case l_s = 0.3 indicates that the original mechanism identified in Ref. 17 is still at play. By visual inspection of this figure, the reader will realize that the flow essentially consists of a single roll, which periodically changes its sense (from the clockwise to the counterclockwise direction and vice versa).

For illustration purposes, Fig. 9 shows separately the flow in the xy and yz plane. It can be easily inferred that the existing main roll has a plane of symmetry inclined by 45° with respect to the yz plane. In line with the mechanisms depicted above,¹⁷ due to the alternating direction of the vibrations, the particles periodically hit the walls perpendicular to the y direction and are accumulated accordingly along these walls. However, due to the presence of the inclined roll, which tends to displace these particles due to viscous (drag) effects, the accumulation surface is bent until two closed surfaces emerge in proximity to each wall. The inclination of the roll and its departure from the classical two-dimensional vortex, which would be produced in a square (2D) cavity (due to the adjustments that its external shape must accommodate in order to fit the cubic symmetry of the considered cavity), are responsible for the double-cone shape of the particle formations.

The next figure of the sequence (Fig. 10), however, shows that this peculiar flow configuration is lost, when the spot size is increased to l_s = 0.4. For simplicity, the particle structures are not shown in this figure (the reader being referred to Fig. 4 for its morphology). In agreement with the theory elaborated by Santhosh and Lappa,²⁵ in this case the size of boundary-dense structures is greatly reduced, while particle-free-boundary formations manifest in the central region. This is due to the existence of stages of evolution [see Figs. 10(b) and 10(d)], where in place of a single roll pervasive throughout the cavity, a two-roll configuration with a smaller extension along y is produced. Accordingly, a relatively long part of the vibration period exists where particles are prevented from reaching the lateral solid walls where they can accumulate and form dense surfaces, which explain the observed decrease in the lateral formations size and the emergence of new free-boundary structures in the center.

A similar situation is obtained in the case l_s = 0.5 (see Fig. 5 for the corresponding particle structures). Notably in this circumstance, the lateral particle-dense-boundary structures are completely suppressed, thereby leaving the floor exclusively to the boundary-free structures. Once again, a detailed assessment of the corresponding fluid-dynamic field is beneficial for the elaboration of a proper explanation. As witnessed by Fig. 11, indeed, in this case no stage of evolution can be recognized where a single roll pervasive throughout the cavity exists. As a result, the mechanism originally identified by Lappa¹⁷ cannot be effective and only particle-free-boundary formations are allowed.

Finally, for l_s > 0.5 (Figs. 6 and 7), as one would expect on the basis of relatively simple arguments, compact boundary-dense particle structures reenter the dynamics because in this range of the spot size, the flow features again a single pervasive roll throughout the duration of the forcing period (not shown).

For the sake of completeness and for consistency with Ref. 25, we also describe the results obtained when the vibrations form an angle ϕ = 45° with the y axis. As the reader will immediately realize by taking a look at the related patterning behavior (shown in Figs. 12–14 for some representative values of l_s), in these cases the formations are relatively irregular, which makes the related treatment somehow worthwhile.

For relatively small values of l_s, as even a cursory comparison of these figures with those reported in Sec. IV A would immediately confirm, their distinguishing mark is a relatively filamentary nature (see, e.g., Fig. 12 for l_s = 0.3).

Moving on to the case with l_s = 0.5 (Fig. 13), it can yet be seen that no well-defined compact structures such as those obtained for ϕ = 0° are visible. Rather, particles tend to accumulate along almost one-dimensional filaments that intertwine and form like an irregular “fabric” or web [Fig. 13(d)].

As a final look at Fig. 14 would indicate, however, well-defined compact (although not symmetric) structures are recovered for l_s ≥ 0.8, which once again might be interpreted as a significant change in the underlying cause-and-effect relationships driving particle clustering (the reader being referred to the final discussions implemented in Sec. IV D).

Following the same approach undertaken in Sec. IV B, in this final section an interpretation of the observed dynamics is attempted on the basis of the possible formation mechanisms elucidated in that section.

Along these lines, we start from the fascinating case shown in Fig. 12 for which some compact “rings” are obtained in place of spatially extended surfaces. Also in this case, direct analysis of the fluid-dynamic behavior qualitatively shown in Fig. 15 is instrumental in getting useful insights into the particle dynamics. On close inspection, this figure reveals that stages of evolution exist where the flow in the xy plane features only one eccentric dominant roll [i.e., located in proximity to one of the corners, see, in particular, Figs. 15(b) and 15(d)] and two other concentrated (small) vortices. Apart from the asymmetry related to the position of this roll, an important aspect distinguishing this case from all the others is the fact that, although the sense of rotation of this roll changes when moving from the time t = t₀ + P/4 to the time t = t₀ + 3P/4, however, its eccentric location remains exactly the same. This might be seen as the sought justification for the emergence of the particle formation just in proximity to this corner. One may assume that the mechanism originally identified by Lappa¹⁷ is effective essentially in this area.

A separate discussion, however, is needed for the formations visible in Fig. 13. We ascribe the involved filamentary nature of the structures and their intertwined appearance to the extremely multicellular nature of the carrier flow, which as illustrated in Fig. 16, is made of different coexisting rolls [up to five at the same time as shown in Fig. 16(c) in the xy plane]. Moreover, in this case, no stage exists at all where a single pervasive roll is established in the cavity.

As a concluding remark for this section, we limit ourselves to highlighting that the compact nature of the structures visible in Fig. 14 stems once again from the ability of the system to produce for this value of the spot size a single roll with given orientation inside the cubic enclosure, which changes periodically its sense of rotation (Fig. 17).

To complement an earlier investigation with new findings naturally stemming from an increase in the multiplicity of the thermal inhomogeneities (in the form of central spots), a new numerical study has been conducted to assess the nature and morphology of the emerging particle structures. These have been found to show up with a kaleidoscope of variants, which depend on the considered size of the spot and the inclination of vibrations.

In particular, the structures come under the heading of two fundamental categories, namely, the compact formations whose distinguishing mark is a well-defined boundary formed by the repetitive collision of particles with the enclosure walls, and alternate aggregates where such a feature is not present as the interaction with the walls does not a play a crucial role.

Both seem to be possible in a certain region of the space of parameters, however, when the vibrations are perfectly perpendicular to the adiabatic cavity walls, the latter type manifests itself only in a very restricted range of spot sizes. As made evident by a detailed analysis of the underpinning velocity field, the morphology of such boundary-free formations depends essentially on the alternate multi-roll configurations, which temporarily replace (during one period of vibrational forcing) the classical single vortex typical of thermovibrational flow. Their size is directly connected to the duration of the timeframe in which such multi-roll convective patterns are effective.

Completely new phenomena are observed in the case for which the vibrations are inclined with respect to the walls. For relatively small values of the spot size (l_s = 0.3 or 0.5), the particle aggregates lose their compact nature and tend to be replaced by more or less filamentary formations, which in some cases are heavily intertwined, thereby giving the observer the illusion of a “fabric.” Such a complexity is supported by evident temporal and spatial symmetry breaking phenomena affecting the topology of the underlying (carrier) fluid flow. In such circumstances, the presence of a single roll, pervasive throughout the cavity, is no longer a flow feature.

To allow for additional impact of this line of research, several steps could be taken in the future to explore even further the consequences of the wall temperature distribution on the system dynamics and consider cavities with shapes more complex than the cubic (fundamental) one.

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 authors have no conflicts to disclose.

Balagopal Manayil Santhosh: Data curation (equal); Software (equal); Visualization (lead); Writing – review & editing (supporting). Marcello Lappa: Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (equal); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are openly available in Pure, at https://doi.org/10.15129/8ad286e1-bf87-42cb-88b9-fc4fe392676e (Ref. 34).