Numerical investigation of induction of chaotic micromixing via vibration switching Open Access

Microfluidic technologies have significantly improved biomedical and chemical reaction technologies for several decades.^1,2 Although mixing is an important operation in chemical reactions, separation, and various biochemical assays,³ fluids are not readily mixed because the microscale flows in microfluidic chips are primarily laminar owing to the low Reynolds number, i.e., the dominance of the viscous force over the inertial force. Therefore, a mixing unit is required in most integrated microfluidic devices, and various passive and active mixing techniques have been developed.^4–6 

The chaotic advection theory has been proposed as a strategy for promoting mixing in a laminar flow regime.^7–9 The essential idea of this scheme is to exponentially increase the interface between layers by stretching and folding the fluid elements (Baker's transformation), even for systems dominated by viscosity. Although only the stretching of fluid elements occurs by shear in steady laminar flows, folding can be implemented by switching the flow field such that the streamlines at different instants intersect. Cavity and blinking vortex flows have been proposed as classical and representative chaotic mixing systems, respectively.^7,10

Inspired by these studies, researchers have attempted to implement chaotic mixing in microfluidic systems. Chaotic mixing cannot be induced in simple steady laminar flows; therefore, various strategies have been investigated. Regarding passive chaotic mixing, for example, microfluidic channels with grooves at an angle slanted toward the flow direction have been designed to induce secondary flows (vortices along the flow axis). The alternate installation of slanted grooves with different angles allowed the rotational direction of these vortices to be switched along the stream, thereby resulting in chaotic mixing.^11,12 Three-dimensional (3D) folding of the channel path has been shown to induce stretching and folding in fluid elements at a moderate Reynolds number.^13,14 Furthermore, efficient mixing achieved through abovementioned techniques has been exploited to develop devices for capturing specific cells (e.g., tumor cells).^15,16 Although these devices have been proven to be effective, they are primarily flow-based systems requiring a pump and complex channel geometry, thus resulting in a bulky system and complex fabrication processes. Meanwhile, active mixers have been designed to promote mixing by applying external inputs such as pressure, electrokinetic, or magnetic forces.^17–19 Active mixing schemes can be applied to both flow- and static-type systems, thus affording a higher degree of flexibility for the mixing methods. Although examples that utilize electricity or magnetism exist, their applicability is limited because of their dependence on the material properties of the fluid or particles, and the systems involved are complex. Therefore, for practical applications, a simple system is required to realize versatile fluid control.

As an active fluid control method, the application of the surface acoustic wave (SAW) to fluid is well recognized, involving the propagation of minute waves on the substrate through resonating electrodes.^20–24 While widely utilized for fluid/particle control without the need for an external pump system, this method requires intricate electrode designs and fabrication procedures. In addition to the SAW technique, steady streaming (SS) is also attracting attention as an alternative. SS can be achieved by introducing sound waves or vibrations from a wave source, enabling localized manipulation of microscale flow fields around microstructures without complex fabrication processes. This approach is straightforward and allows for flow control via various parameters, such as the vibration frequency, amplitude, and geometries of microstructures.^25–31 In the context of SS-based micromixers, Huang and co-workers^32–34 demonstrated mixing based on a localized flow induced around a sharp-edge. Additionally, chaotic mixing within microwells using audio speakers has been demonstrated, even in cases where microstructures are absent.³⁵ 

Although SS is characterized by its simplicity in terms of system and operation, predicting its behavior is not straightforward. Thus far, two-dimensional (2D) theoretical analyses and numerical simulations have been primarily conducted.^36–39 Additionally, 2D numerical modeling was performed for the design optimization of the abovementioned sharp-edge micromixers to investigate the effects of multiple parameters, such as the sharp-edge angle and vibration amplitude.^40–43 However, because most microfluidic devices are shallow in the depth direction, the quasi-2D approximation is not valid and consideration of the 3D velocity distribution is necessary. Although experimental approaches can clarify the characteristics of 3D flows, they require expensive equipment, such as confocal microscopes, to conduct fine-resolution 3D particle tracking velocimetry (PTV). Therefore, the ability to numerically reproduce 3D flows is key for achieving a comprehensive understanding of SS flows.

Thus far, we have investigated the application of vibration-induced flow (VIF) to micromixing. To obtain the design guidelines for the mixer, we developed a 3D numerical simulation tool,^28,31 which allowed us to predict various flow fields for the flexible adjustment of parameters. In this study, we hypothesize that flow field switching implemented by periodically altering the vibration condition can result in chaotic advection. We verify this idea by examining the flow field and trajectories of fluid particles using a previously developed 3D numerical simulation tool. Additionally, we quantitatively assess the extent of chaotic advection using the Lyapunov exponent.^44,45

In VIF, four localized vortices are generated around a pillar when a rectilinear vibration is applied in the horizontal direction.^27,36–38 The positions of the vortices are determined based on the vibration axis. In this study, we simulate the flow field induced by rectilinear vibrations applied to a substrate with micropillars in vibration directions of 0°, 45°, 90°, and 135° with respect to the x axis (Fig. 1). Subsequently, we assess whether the switching of the vibration directions distorts the fluid element around the pillar. By combining these different flow fields in periodic switching, we expect to realize chaotic advection based on Baker's transformation via simple vibrations.

The flow fields around a micropillar placed between two parallel plates are calculated by assuming that the liquid is a Newtonian fluid (water) and incompressible. The computational domain is shown in Fig. 2, where the top and bottom boundaries correspond to two parallel plates. A pillar is fixed to the two plates at the center of the computational domain. The half-height of the computational domain δ is set to 50 μm. The two orthogonal directions tangential to the substrate are set as x and z, whereas the wall-normal direction is y. The origin is located at the center of the x–z plane. We also examined the effect of the pillar shape by comparing cylindrical and triangular prism pillars. For the cylindrical pillar, the diameter and height of the pillars are set to be 100 μm (2δ), which is the same as the height of the computational domain [Fig. 2(a)]. For the triangular pillar, the length of one side is set to 100 μm, and the height is the same as the cylindrical pillar [Fig. 2(b)]. In the computational domain, no-slip conditions are imposed at the bottom and top plates. Periodic boundary conditions are imposed in the horizontal x- and z-directions, indicating that the pillars exist in a regular matrix (array). The effects of pillar geometry and vibration conditions on chaotic advection are investigated by varying the pillar-to-pillar separation S and vibration amplitude A. Specifically, we set the above parameters in the range of S = 150 (3δ), 200 (4δ), and 400 μm (8δ), and A = 1.5, 3.3, and 4.5 μm, with a constant frequency f = 600 Hz.

In our simulation, a moving coordinate system that propagates with the plates and pillar is considered instead of a stationary coordinate system to avoid the calculation of moving boundaries. The effects of vibrations are incorporated by applying a time-varying horizontal body force (inertial force) to the fluid region. The flow field obtained in the moving coordinate system is then converted to one based on a fixed coordinate system.

We calculated the flow fields induced by rectilinear vibrations in four different directions [0°, 45°, 90°, or 135° with respect to the x axis; Fig. 2(c)]. The formulas expressing the displacements of the solid substrate for each rectilinear vibration are shown in Fig. 2(c), where t and ω represent the time and angular frequency of vibration ω = 2π/T, respectively (T represents the vibration period).

In the preliminary investigation, we found that the induced velocity magnitudes for the 45° and 135° axes are smaller than those for the 0° and 90° axes when an identical vibration amplitude A was employed (Fig. S1 in the supplementary material). Thus, we employed the composite amplitude $A ′=A× 2$ instead of A for the 45° and 135° axes to match the induced velocities for each vibration axis in the vibration-switching simulation.

In our calculations, the equations above are solved after normalizing all physical quantities with δ and T.

The governing equations shown in Eqs. (1) and (2) in the moving coordinate are solved using a pseudo-spectral method, in which a solution is expanded by Fourier modes in the x- and z-directions, and Chebyshev polynomials in the y-direction. For time advancement, the Crank–Nicolson scheme is used for the diffusion terms, whereas the second-order Adams–Bashforth scheme is used for the convection terms.^46,47 The number of modes employed in the current simulation are (N₁, N₂, N₃) = (64 × 33 × 64) in the x-, y-, and z-directions, respectively. The 3/2 rule is used to remove aliasing errors such that the nonlinear terms can be evaluated at 1.5 times finer physical grid points in each direction. The numerical time step, which is nondimensionalized by T, is set as Δt* = 1.0 × 10⁻⁴. The superscript asterisk indicates a dimensionless quantity.

In our simulation, we use a volume penalization method to apply a no-slip condition to the solid surface of the pillar, which enables a complex structure to be embedded in the Cartesian coordinate system. The geometry of the pillar is first expressed by a level-set function ϕ₀,⁴⁸ which is then converted to the phase-identification function ϕ. ϕ is set to be unity in the solid and null areas of the fluid regions. The phase-identification function changes seamlessly from 0 to 1 across the interface within a few grid points, thus avoiding numerical instability. The no-slip condition at the solid surface is imposed by introducing an artificial damping force to the Navier–Stokes equation proportional to ϕ. Calculation is performed until the velocity field induced by a periodic inertial force is sufficiently converged. Finally, 50 instantaneous Lagrangian velocity fields around the pillar are obtained under each condition.

In VIFs, a fluid particle merely oscillates and no net flow arises in the void space; however, the presence of solid microstructures in the domain slightly shifts the position of the fluid particle when it returns to the same phase during periodic oscillation, thus resulting in a non-zero net flow. To capture this phenomenon, a time-averaged Lagrangian velocity field is calculated to characterize the features of the flow field. We performed Lagrangian tracking of convecting fluid particles based on a periodic instantaneous Lagrangian velocity to obtain a net flow field in the fixed coordinate system. Similar to our previous studies,^28,31 we calculate the trajectories of fluid particles initially distributed at regular grids in a 3D computational domain using the fourth-order Runge–Kutta method¹⁹ with a time step of Δt* = 5.0 × 10⁻³. In this process, the instantaneous velocity between the grids is linearly interpolated from the velocities of the eight surrounding grids. Additionally, linear interpolation in the time dimension is conducted using the velocity fields of two neighboring frames. After tracking the fluid particles for ten periods of vibration, the time-averaged Lagrangian velocity field is obtained from the displacement between the start and end points of each particle.

To examine the characteristics of vibration-induced chaotic advection, a Lagrangian fluid particle tracking simulation is conducted using an instantaneous (Lagrangian) flow field. The time-averaged Lagrangian velocity presented in Sec. II D represents the particle motion with time averaging. However, in reality, fluid particles propagate with periodic oscillations at the instantaneous velocity (an example of the actual trajectory of a particle and the derivation of the time-averaged Lagrangian velocity is shown in Fig. S2 in the supplementary material). Therefore, to precisely reproduce the behavior of the advection path of a fluid particle, instantaneous velocity fields must be used instead of time-averaged Lagrangian velocity fields. In this study, we performed particle tracking for particles initially distributed in specific regions in the 3D fluid domain using the same numerical method used for Lagrangian time averaging. To visualize the particle trajectories, the particle positions at identical phases in the time series of periodic vibrations are exported and plotted.

The extent of chaotic advection is evaluated using two approaches based on fluid particle tracking. First, the deformation behavior of fluid particle lumps is visualized. Two lumps of fluid particles (with each lump comprising 6400 particles) are located near the pillar with a vertical position of y = 47.5 μm (the actual initial position of the lumps is shown in Fig. S3 in the supplementary material). The particle spacing in the x–z direction was set to 0.25 μm.

First, we examined the time-averaged Lagrangian velocity field induced by rectilinear vibration around the cylindrical pillar. Figures 3(a)–3(d) show 2D vector plots of the calculated velocity fields at the middle plane of the computational domain (y = 47.5 μm) for four vibration axes. Under rectilinear vibration in the x axis [0° vibration; Fig. 3(a)], four vortices are formed in each quadrant around the pillar, with the center located at a 45° offset from the vibration axis, which is similar to the findings of previous studies.^27,36–38 Consequently, the inflow and ejected flows are induced parallel and perpendicular to the vibration directions, respectively. The velocity field induced by the vibration along the z direction [90° vibration; Fig. 3(c)] is 90° rotationally identical to that along the x-direction, as expected. The velocity magnitude in the fast region is approximately 100 μm/s in each vibration direction, and the velocity of the ejected flow exceeds that of the inflow. Under 45° and 135° vibrations [Figs. 3(b) and 3(d), respectively], the vortices are diagonally offset from those under the 0° and 90° conditions. In addition, for the 45° and 135° vibrations, the velocity of the inflow exceeds that of the ejected flow, in contrast to the cases for 0° and 90° vibrations. We speculate that a higher inflow velocity is induced because the inflow path toward the pillar wall is narrower owing to the shift in the vortex positions.

Figures 3(e) and 3(f) show 2D vector plots in the vertical (y) sections under 0° and 45° vibrations. The top and bottom panels show planes parallel (x and $C− C ′$ sections) and perpendicular (z and $D− D ′$ sections) to each vibration direction, respectively. Under 0° vibration, inflow toward the pillar wall is induced in the plane-parallel to the direction of the vibration [Fig. 3(e); top panel]. By contrast, an ejected flow emerges in the plane perpendicular to the vibration direction [Fig. 3(e); bottom panel] with a higher velocity than that of the inflow. Under 45° vibration [Fig. 3(f)], both inflow and ejected flows are induced in planes parallel and perpendicular to the vibration direction. However, in contrast to the 0° condition, the inflow exhibits a higher velocity compared with the ejected flow.

In previous studies involving 2D analyses,^27,36 the velocity along the vortices circulating in two dimensions is expected to remain relatively consistent. However, the 3D simulation results calculated using our current code show differences between the inflow and ejection flows. This indicates that the flow field around the pillar is affected by the 3D behavior of the flow.

Figure 4 shows representative snapshots of two lumps of fluid particles (denoted in red and blue; actual movies are shown in Movie S1 in the supplementary material). Differences in the sequence of vibration switching are compared to determine the mechanism by which the switching of the flow field affects the deformation of the fluid elements at a switching interval T_sw = 0.83. At this T_sw, St is approximately 1; therefore, the advection of the fluid in this period is comparable to the scale of the flow field. In the case of no switching as a reference [Fig. 4(a)], only stretching deformation of lumps of fluid particles along the local vortices is induced (6.3 and 12.5 s). Under the switching of the 0° and 90° axes [Fig. 4(b)], the lumps stretch and fold (8.8 s), and the interface between the layers of the lumps increases. Nevertheless, the particles remain within the small regions and do not propagate to the wide flow domain (12.5 s). This should be because only the rotational direction of the vortices is reversed by switching the flow fields [Figs. 3(a) and 3(c)]. Under the switching of the 0° and 45° axes [Fig. 4(c)], the lumps stretch along the flow first (8.2 s), and once the flow field is switched, the particles in a lump are drawn into two different vortices and segregated into two elements (9 s). By contrast, under four-axis switching [Fig. 4(d)], the lumps are stretched following a flow toward the pillar wall (6.6 s). Subsequently, as the flow switches, the particles are drawn into different vortices and stretched (7.4 s). In the latter two cases [Figs. 4(c) and 4(d)], the lumps of particles stretch and fold; simultaneously, they are advected in the entire computational domain. These results indicate that periodically switching different flow patterns and changing the vortex positions concurrently can effectively distort the lumps of fluid particles and transport them over a wide area.

Next, we calculated the long-time trajectories of single particles for 200 s to determine the spatial range in which single elements of the fluid can be examined under conditions identical to those shown in Fig. 4. Figure 5 shows a trajectory projected in the x–z plane and the time course of the y position (t–y plot) of a representative single particle that initiates from (x, y, z) = (59.9, 47.5, 59.9). In the no-switching case [Fig. 5(a)], the particles propagate along the vortex flow in a quadrant, which is similar to the lump deformation shown in Fig. 4(a). For the t–y plot in the bottom panel, the particle oscillates quasi-periodically at y = 50–100 μm. Under the two-axis switching of 0° and 90° vibrations [Fig. 5(b)], a particle repeatedly propagates back and forth along a path of vortices in reverse directions as the flow fields switch, and the particle remains within a small exploration range. Similar to the no-switching case, a particle shows quasi-periodic 3D motion, thus indicating that the particle propagates not only in the horizontal (projection to x–z plane) direction but also in the vertical direction [bottom panel of Fig. 5(b)]. In both the 0° and 45° axes and the four-axis switching conditions [Figs. 5(c) and 5(d)], a particle explored the entire x–z plane by transitioning four quadrants. The stretching and folding of the lumps shown in Figs. 4(c) and 4(d) are caused by this active particle advection. Additionally, the particles exhibit irregular motion in the vertical direction under both two-axis of 0° and 45° vibration and four-axis conditions. These results suggest that the periodic switching of the vibration directions contributes to chaotic motion in both the horizontal and vertical directions.

For the aforementioned vibration conditions, the LE was calculated to evaluate the extent of chaotic advection. Figure 6 shows the time evolution of the LE values for each vibration condition. In the cases of no switching (black line) and two-axis of 0° and 90° vibrations (blue line), the LE converges to approximately zero after 200 s. Meanwhile, it converges to approximately 0.4 under the two-axis of 0° and 45° (yellow line) as well as four-axis vibrations (red line). The behavior of the LE value reflects the visualization results of lump deformation and single fluid particle tracking. High LE values are indicated in cases with intense advection and better mixing in the computational domain, as seen in Figs. 4 and 5. This correlation between the LE and lump deformation, which was similarly observed in previous chaotic mixing systems, is confirmed in our system. Based on the resulting particle trajectories and LE values, we can conclude that by employing the switching of vibrational directions, features of chaotic advection and mixing emerged even when simple cylindrical pillars are adopted.

Next, we investigated the dependence of the vibration-switching interval (i.e., the effect of St) on chaotic advection with two-axis of 0° and 45° vibrations. Figures 7(a)–7(c) show images of the lumps of fluid particles at 12.5 s, with the initial locations identical to those presented in Fig. 4 (actual movies are shown in Movie S2 in the supplementary material). Under a short switching interval [T_sw = 0.08 s (St = 9.5), Fig. 7(a)], the lumps of the fluid particles are perturbed only within a small area. Under a long switching interval [T_sw = 8.3 s (St = 0.095), Fig. 7(c)], the stretching period is too lengthy, which hinders the efficient induction of fluid lumps movement across the entire computational area. Meanwhile, under an intermediate switching interval [T_sw = 0.83 s (St = 0.95), Fig. 7(b)], the fluid lumps deform in a complex manner owing to repeated stretching and folding, as seen in Fig. 4(c).

Figure 7(d) shows the relationship between St and LE obtained at 200 s for the different switching patterns. Under the switching of 0° and 90° vibrations (blue line), the LE values are as low as 0–0.1 within a St range of 0.05–5. By contrast, under two-axis of 0° and 45° vibrations (yellow line) and four-axis vibration (red line), the LE varies depending on changes in St, with a maximum LE of approximately 0.4. recorded at St = 0.5. These results quantitatively clarify that the extent of chaotic advection depends on the St defined by the switching interval of the vibrations. In the following calculations, we employ the two-axis of 0° and 45° switching condition, which yields relatively high LE values via a simple switching procedure.

Next, the effects of various vibration conditions on the chaotic advection were investigated. Specifically, we examined the effect of the vibration amplitude, separation distance between pillars, and cross-sectional shape of the pillar under the condition of vibration switching shown in Figs. 4(c) and 5(c) as a reference. Before examining the behavior of fluid deformation and the LE, we examined the time-averaged Lagrangian velocity field in the middle (y = 47.5 μm) plane (Fig. 8). The size of each panel in the figure is determined by the corresponding pillar separation, as the figure depicts conditions under various pillar separations. As the vibration amplitude increases [A = 4.5 μm, Fig. 8(a)], the magnitude of the induced velocity almost doubles, whereas the overall shape of the flow field almost resembles that of the reference case [A = 3.3 μm, Figs. 3(a) and 3(b)]. Four vortices are formed in four quadrants with respect to the axis of the vibration direction, as similarly shown in Figs. 3(a) and 3(b). Next, we examined the effects of smaller and larger separation distances between pillars of identical sizes. Under a smaller separation distance [S = 150 μm; Fig. 8(b)], the magnitude of the induced velocity almost halves (∼50 μm/s). We speculate that this is because the interference between vortices induced by adjacent pillars becomes stronger and thus suppresses the velocity magnitude (note that a periodic boundary condition is introduced in the x- and z-directions). Meanwhile, under a larger separation distance [S = 400 μm; Fig. 8(c)], the size of the vortices is smaller than that of each quadrant in the computational domain, which implies that flows at adjacent pillars do not interact. At a larger pillar separation distance (S = 600 μm; Fig. S8 in the supplementary material), the vortices induced around the pillar behave independently, with almost no interactions between those induced by adjacent pillars. Under the 45° vibration condition [right panel of Figs. 8(a)–8(c)] in the aforementioned pillar separation distances, the positions of the vortices are diagonally offset from those at 0°. Next, we evaluated the effect of the noncircular sectional shape using a triangular prism pillar with A = 3.3 μm and S = 200 μm [Fig. 8(d)]. High velocity (∼250 μm/s) flows ejecting from the sharp-edges in a direction perpendicular to the vibrational direction were observed. Accordingly, a pair of vortices are generated around the ejected flow. In the case of 45° vibration [right panel of Fig. 8(d)], the occurrence of high-velocity ejected flows originating from the sharp-edges at two locations increases. Along with the positional change in the ejected flow, the positions of the vortices shift.

Figure 8(e) shows the relationship between the peak magnitudes of the time-averaged Lagrangian velocity fields under the conditions examined in this section [the hatched bars in each panel represent the reference case shown in Figs. 3(a) and 3(b)]. The magnitude of the induced velocity increases with the vibration amplitude (left panel). Similarly, the velocity increases with the pillar separation distance (center panel). The flow field induced around the triangular prism pillar shows a higher velocity that is thrice that of the cylindrical pillar under identical vibration conditions (right panel) because of the ejection flow emitted from the concentrated area around the edge of the pillar, as reported previously.^32,40–42

The extent of the chaotic advection was quantified by evaluating the LE based on the obtained flow fields under various vibration conditions. Here as well, an instantaneous flow field is used to track the particles in the flow induced by a two-axis switching (0° and 45°) condition. Figures 9(a)–9(e) show the trajectories of single tracer particles for 200 s projected onto the x–z plane of each geometry. Under all conditions, a particle traverses a large portion of the flow domain. Comparing conditions with different amplitudes [Figs. 9(a) and 9(b)], a particle advected more actively in the higher amplitude condition [A = 4.5 μm, Fig. 9(b)]. Additionally, in both cases with smaller [S = 150 μm, Fig. 9(c)] and larger [S = 400 μm, Fig. 9(d)] pillar separation distances, a particle is advected throughout the entire x–z plane. In the case of the triangular prism pillar [Fig. 9(e)], a particle propagates actively in the entire x–z plane, which is similar to the case of the cylindrical pillar with a larger vibration amplitude [A = 4.5 μm, Fig. 9(b)]. We observed that the 3D characteristics of the trajectories differ among the conditions. The t–y plots of the tracer particle under each condition (Fig. S9 in the supplementary material) clearly elucidate these behaviors. In the reference case [Fig. S9(a) in the supplementary material], a particle initially remains within the middle plane but abruptly fluctuates significantly in the vertical (y) direction after 100 s. Under a larger amplitude [Fig. S9(b) in the supplementary material], a particle oscillates at an irregular frequency. The y position of a particle fluctuates significantly under the case involving a smaller pillar separation distance [S = 150 μm, Fig. 9(c)]. However, under a larger pillar separation distance [S = 400 μm, Fig. 9(d)], the particle remained within the middle plane. Based on these results, the emergence of irregular trajectories in both the horizontal and vertical directions implies the occurrence of chaotic advection.

Figure 9(f) shows the relationship between the LE values after 200 s for the abovementioned vibration conditions. Regarding the amplitude dependence (left panel), the LE value increases with the amplitude. Among the cases of different pillar separation distances (center panel), a relatively high LE of ∼0.4 is achieved under S = 200 and 400 μm, whereas LE values of 0.2 and 0.3 are achieved under S = 150 and 600 μm, respectively. These results indicate that the LE values increase with the velocity magnitudes shown in Fig. 8(e), except that the LE value decreases under the largest pillar separation distance (S = 600 μm). We speculate that this is because the particle enters a low-velocity region far from the pillar under larger pillar separation, wherein the particle motion decelerates. For the case with the smallest pillar separation distance (S = 150 μm), the interference of vortices becomes apparent because the separation distance between the pillars is too small, which consequently suppresses the particle motion due to the lower velocity. Therefore, the appropriate separation distance between pillars should be selected to achieve higher LE values. Furthermore, the LE for the case with the triangular prism pillar is approximately 1.5 times greater than that for the case with the cylindrical pillar under the identical vibration condition (right panel). This might be because the particle advection is promoted by the localized high-velocity flow near the pillar edges. Finally, the St–LE plots for all evaluated conditions and geometries (Fig. S10 in the supplementary material) show peaks at St ∼ 1. This indicates that optimal chaotic advection can be achieved when the timescales of advection and flow switching are comparable to induce an efficient exponential lamination of fluid elements.

Finally, the 3D behavior of the fluid particle lumps is examined via particle tracking. A total of 7200 particles are initially placed in the horizontal planes, with y = 47.5 and 45 μm (actual initial positions are shown in Fig. S11 in the supplementary material). Figure 10 shows snapshots of the particle tracking simulations at 12.5 s for each vibration condition (actual movies are shown in Movie S3 in the supplementary material). The top and bottom panels represent projections onto the x–z and x–y planes, respectively. Figures 10(a) and 10(b) depict the cases for A = 3.3 and 4.5 μm, respectively, showing that the initially localized particles are distributed more uniformly in the fluid domain in the latter case. Based on the side views (bottom panels), a higher amplitude promoted more intensive 3D advection of the particles within the lower half of the fluid domain. Figures 10(c) and 10(d) show the results under S = 150 and 400 μm, respectively. Under S = 150 μm, the particle lumps remain concentrated within the small region. Under S = 400 μm, particles traverse farther from the pillar with almost no 3D motion. For the case with the triangular prism pillar [see Fig. 10(e)], the particles are advected similarly to the cylindrical pillar with a larger vibration amplitude [A = 4.5 μm, Fig. 10(b)].

Based on the results shown in Fig. 10, two types of particle lumps are well mixed, particularly under the case involving A = 4.5 μm [Fig. 10(b)] and the triangular prism pillar [Fig. 10(e)]. These two conditions are associated with high LE values, as confirmed based on Fig. 9(f); this indicates that the LE reflects the mixing performance, which agrees well with the findings of previous studies.^17,19,45 The active particle advection observed under these conditions is caused by Baker's transformation of the stretching and folding of the fluid element (Fig. 4). Therefore, we conclude that switching the VIF is effective in inducing efficient mixing.

We investigated the effect of periodically switching VIF fields on the induction of chaotic advection by combining simple rectilinear vibrations in different directions. The novelty of this study lies in its extension of our previously reported numerical model,^28,31 enabling us to elucidate the behavior of fluid particles and the deformation dynamics of particle lumps via a particle tracking simulation based on a periodic instantaneous flow field in 3D space; we performed such a simulation as experimental observations would not be straightforward. The simulation allows us to determine the mixing characteristics that depend on geometric and vibration parameters.

Additionally, achieving chaotic mixing in a simple system by periodically switching the flow induced by rectilinear vibrations is also a novelty. By periodically switching the flow fields, we have observed the stretching and folding of lumps of fluid particles, which resulted in temporally and spatially irregular advection, as observed in classical models of chaotic mixing, such as those for cavity flow¹⁰ and blinking vortex flow.⁷ Moreover, our 3D simulation results show that intensive 3D advection might occur in VIFs with quasi-2D geometries depending on the conditions, whereas most previous studies discussed only 2D regular (steady) behaviors.^27,36–39 Our present results indicate the possibility of 3D complex mixing patterns arising from simple geometries and driving modes.

We quantified the extent of chaotic advection using the LE. The irregular chaotic behavior of fluid particles, which was characterized by the stretching and folding of fluid lumps, resulted in a higher LE value, thus confirming the previously reported relationship between the LE and chaotic behavior. Furthermore, the LE values depend on St, which characterizes the ratio of the time scales of the mean flow and the vibration-switching period. Therefore, chaotic mixing can be enhanced by selecting an appropriate switching period.

Additionally, we examined the effects of vibration conditions such as the vibration amplitude, pillar separation, and cross section of the pillar. The general conclusion obtained is greater magnitudes of the induced flow velocities result in a better mixing state within a certain period. However, this principle is not applicable to larger pillar separation distance. In the VIF, the region of nonzero net flow is confined around a micropillar within the spatial range characterized by the Stokes layer thickness. Therefore, the low-velocity region enlarges under larger interpillar separation distances, and the overall mixing performance in the entire space deteriorates. Meanwhile, under too small separation distances, the magnitudes of the induced velocity decrease owing to the interference of adjacent vortices, which consequently suppress the fluid particle motion. Hence, an optimal pillar separation distance that enhances particle exploration and increases the LE value exists. Furthermore, the micropillar with sharp-edges generates strong localized flows, thereby enhancing particle advection, unlike a micropillar with a round edge. The LE of our mixing system reaches a maximum value of approximately 0.8, which is comparable to or higher than the values reported in previous studies that implemented chaotic advection in microfluidics.^17,19,45 Visualization of fluid particles in 3D space revealed that particle lumps traverse three-dimensionally. Although this 3D mixing resulted from the simple 2D vibrations could be a promising approach to develop an efficient micromixer, the advection is limited within either the upper or lower half of the computational domain. We speculate that an extensive mixing across the entire 3D space could be achieved by remains introducing micropillars with 3D-shapes, which will remain for our future studies.

In summary, we numerically demonstrated the occurrence of chaotic advection and the resultant mixing enhancement based on the VIF. We showed that even a simple cylindrical pillar can generate chaotic behavior via the periodic switching of the rectilinear vibration directions. Our numerical tool, which enables us to predict the behavior of fluid particles in 3D space, is useful for optimizing mixing facilitated by localized flow around oscillating microstructures. Especially, the primary application envisioned is the mixing of macromolecules and nano/micro particles with low diffusion coefficients. Furthermore, because VIF techniques have been shown invaluable for other fluidic manipulations, such as liquid/particle transport,^25,50 sorting,^51,52 trapping,^27,36–39 and alignment,^26,53 we plan to combine various units based on VIF control to realize a pump-free integrated lab-on-a-chip system.

See the supplementary material, which includes figures and videos, for additional details regarding flow fields and LE calculation based on particle tracking. Movie S1, Movie S2, and Movie S3 depict actual motions of lumps of fluid particles described in Figs. 4, 7, and 10, respectively.

This study was supported by JSPS KAKENHI Grants (Nos. 19H02115, 19H02576, 20H05935, 19H00901, 22H01454, and 23KJ1940). This study was also supported by the Institute of Science and Engineering at Chuo University, Japan. We would like to thank Zhitai Huang and Yuto Asada for their assistance with flow field analysis and particle tracking.

The authors have no conflicts to disclose.

Kanji Kaneko: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Yosuke Hasegawa: Methodology (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). Takeshi Hayakawa: Funding acquisition (equal); Supervision (equal); Writing – review & editing (equal). Hiroaki Suzuki: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.