Numerical and Experimental Analyses of Three-Dimensional Unsteady Flow around a Micro-Pillar Subjected to Rotational Vibration

The steady streaming (SS) phenomenon is gaining increased attention in the microfluidics community, because it can generate net mass flow from zero-mean vibration. We developed numerical simulation and experimental measurement tools to analyze this vibration-induced flow, which has been challenging due to its unsteady nature. The validity of these analysis methods is confirmed by comparing the three-dimensional (3D) flow field and the resulting particle trajectories induced around a cylindrical micro-pillar under circular vibration. In the numerical modeling, we directly solved the flow in the Lagrangian frame so that the substrate with a micro-pillar becomes stationary, and the results were converted to a stationary Eulerian frame to compare with the experimental results. The present approach enables us to avoid the introduction of a moving boundary or infinitesimal perturbation approximation. The flow field obtained by the micron-resolution particle image velocimetry (micro-PIV) measurement supported the three-dimensionality observed in the numerical results, which could be important for controlling the mass transport and manipulating particulate objects in microfluidic systems.


Introduction
The hydrodynamic phenomenon known as steady streaming (SS) is gaining increased attention for controlling flows and associated transport and mixing of chemical species as well as micro objects such as functionalized particles and cells in microfluidic devices [1][2][3][4][5][6][7][8][9][10][11][12]. This term represents the time-averaged non-zero mean flow induced by relative periodic oscillation with zero mean between the substrate and the adjacent bulk fluid. When an infinite planar substrate is oscillating in parallel to its surface, only a transient velocity field with zero mean is generated within a thin layer whose length scale is characterized by the Stokes boundary layer thickness represented as δ s~( 2ν/ω) 1/2 , where ν and ω are the kinematic viscosity and angular frequency, respectively [13]. However, when an obstacle is present in the flow field, the interaction between the oscillating bulk fluid and the obstacle creates vorticity and causes net-momentum transfer. This is similar to the Reynolds shear stress that arises from the correlation of the fluctuating velocity components in turbulent flows, and it appears as an additional forcing term in the averaged Navier-Stokes equations. As a result, the steady time-averaged velocity is generated, although the applied periodic forcing does not have a mean component. Because it requires no net displacement or a pressure gradient to drive the flow, the SS is expected to simplify and miniaturize microfluidic systems without introducing external pumps or tubing. the boundary is moving relative to the stationary reference frame. The governing equations of an incompressible fluid, i.e., Navier-Stokes and continuity equations, in which a temporally periodic inertia force due to the circular vibration was included, were directly solved by a pseudo-spectral method [24]. Once the numerical results in the moving (Lagrangian) coordinate were obtained, the flow field was converted to the stationary (Eulerian) coordinate. This approach enabled us to calculate the SS flow field without imposing the moving boundary nor the small perturbation approximation. A fluid-solid boundary is expressed by the level-set function [25], which is defined as a signed distance function from the surface. This allows to immerse an arbitrary 3D shape in the Eulerian coordinate system and a no-slip condition at a solid surface is achieved by a volume penalization method (VPM) [26]. Such an approach has a great advantage in implementing arbitrary complex obstacles in the fluid domain without generating numerical grids for each geometry. After obtaining the periodically varying instantaneous velocity field, the time-averaged velocity field of Stokes drift was obtained by tracking the virtual fluid particle imposed within the fluid.
The 3D SS velocity field obtained by the numerical simulation was validated by quantitatively comparing with the 3D micro PIV measurement results obtained from the confocal microscopy equipped with a high-speed camera. Our results exhibited the good agreement even within the steeply varying velocity profile within the Stokes layer. The 3D paths of tracer particles induced by the vortex structure adjacent to the edge of the cylinder observed in the experiment were reproduced in the numerical simulation.

Computational Domain and Governing Equations
In this study, we calculated the flow field around a cylindrical micro-pillar placed between two parallel plates (Figure 1a). Both plates and a pillar fixed to the bottom plate periodically oscillate following a circular path parallel to the substrate (Figure 1b). To simulate this system, we employed a moving coordinate system that moves with the plates instead of the stationary coordinate system with moving boundaries. The flow field obtained in the moving coordinate system was eventually converted to the coordinate system at rest. The advantage of this approach is that the calculation code is simple and numerical accuracy is high, since boundaries (the substrate and the pillar) are stationary with respect to the coordinate system. The effect of the circular vibration was considered by applying an inertia force that rotates in accordance with the acceleration/deceleration of the moving coordinate.
Assuming that the liquid is incompressible and a Newtonian fluid, the liquid flow is governed by the following Navier-Stokes and the continuity equations: Here, variables with an asterisk represents a dimensional quantity. The fluid kinetic viscosity and density are denoted by ν and ρ, respectively. The two orthogonal directions tangential to the substrate are set as x and z, while the wall-normal direction is y. The origin is located at the center of the bottom substrate. Time is defined as t. The under-bar indicates a physical quantity in the coordinate system at rest, which is introduced so as to be distinguished from the moving coordinate along a circular orbit defined later. The fluid velocity and static pressure are denoted by u i and p, respectively, where the subscript of i represents the three directions, namely, i = 1, 2, 3 corresponds to x, y, z, respectively. A micro-pillar was attached to the bottom substrate and its height was set to be equal to the channel half height δ * . No-slip conditions were applied at the bottom and top walls, i.e., y * = 0 and 2δ * , as well as the surface of the pillar. Periodic boundary conditions were employed in the x and z directions.
This condition corresponds to the case where the geometry shown in Figure 1a repeats in these two directions. The present configuration is chosen because it is relatively simple (the two parallel walls and the periodic placement of cylinders between the walls), while the truncation of the cylinder at the middle of the channel causes complex three-dimensionality of the resulting flow. Assuming that the liquid is incompressible and a Newtonian fluid, the liquid flow is governed by the following Navier-Stokes and the continuity equations: * * + * * * = − 1 * * * + * * * * (1) * * = 0 Here, variables with an asterisk represents a dimensional quantity. The fluid kinetic viscosity and density are denoted by ν and ρ, respectively. The two orthogonal directions tangential to the substrate are set as x and z, while the wall-normal direction is y. The origin is located at the center of the bottom substrate. Time is defined as t. The under-bar indicates a physical quantity in the coordinate system at rest, which is introduced so as to be distinguished from the moving coordinate along a circular orbit defined later. The fluid velocity and static pressure are denoted by and p, respectively, where the subscript of i represents the three directions, namely, i = 1, 2, 3 corresponds to x, y, z, respectively. A micro-pillar was attached to the bottom substrate and its height was set to be equal to the channel half height δ * . No-slip conditions were applied at the bottom and top walls, i.e., y * = 0 and 2δ * , as well as the surface of the pillar. Periodic boundary conditions were employed in the x and z directions. This condition corresponds to the case where the geometry shown in Figure 1a repeats in these two directions. The present configuration is chosen because it is relatively simple (the two parallel walls and the periodic placement of cylinders between the walls), while the truncation of the cylinder at the middle of the channel causes complex three-dimensionality of the resulting flow.
In the present study, the solid substrate with micro-pillars oscillates along a given circular tangential orbit as shown in Figure 1b. When the radius and period of the circular oscillation are expressed by A * and T * , the relative displacement of the substrate x * R is given as In the present study, the solid substrate with micro-pillars oscillates along a given circular tangential orbit as shown in Figure 1b. When the radius and period of the circular oscillation are expressed by A * and T * , the relative displacement of the substrate x * R is given as Accordingly, the velocity of the solid substrate, u * R , can be obtained from the time derivative of x * R as Solving Equations (1) and (2) under the boundary condition of Equation (4) on the surfaces of the top and bottom substrates and the pillar requires the treatment of moving boundaries. In order to avoid it, a new coordinate system x * = (x * , y * , z * ) T , which moves with the same speed as the substrate, is introduced. With the reference length scale of δ * and the oscillation period of T * , the generalized dimensionless forms of Equations (1) and (2) on a rotating frame are expressed as where the fluid velocity relative to the substrate motion is defined as u i = u i − u R i . The effect of the rotational vibration appears as a dimensionless inertia force f i due to the acceleration/deceleration of the moving coordinate, which is given by The detailed derivation of Equations (5)-(7) can be found in the Supplementary Materials. Here, the Strouhal number St is defined as the ratio of the time scales of the flow and the oscillation. Considering that the time scale of flow is given by δ * /U * max , where U * max = 2πA * /T * is the maximum velocity of the substrate, the Strouhal number is given by The Reynolds number is given by The above two dimensionless parameters characterize the flow field considered in the present study. Since we solved the velocity field in the reference frame moved with the substrate, the top and bottom walls as well as the micro-pillar stay at rest, and therefore all the velocity components on these boundaries became null. In order to impose a no-slip condition at a solid surface with arbitrary geometry, we used a volume penalization method introduced in the next subsection.

Volume Penalization Method
One of the main objectives in the present study is to develop a numerical code that is capable of simulating a flow around pillars with arbitrary shapes. A volume penalization method is a kind of immersed boundary techniques, in which a complex structure is embedded in the Cartesian coordinate system. In contrast to using a boundary-fitted coordinate system, the immersed boundary technique has the advantages that grid generation is quite straightforward regardless of the complexity of the geometry and highly accurate discretization schemes developed for the Cartesian grid system can be applicable.
In the present study, the geometry of the pillar was first expressed by a level-set function in the Cartesian coordinate system. The level-set function φ 0 is a signed distance function from a surface [25] and it has been widely used to represent complex interface geometry. Here, φ 0 was defined positive inside the solid, and negative in the fluid domain. Then, we converted φ 0 to the phase-identification function φ, which was φ = 0 inside the fluid, whereas φ = 1 in the solid. In order to avoid numerical instability, the phase-identification function smoothly changes from zero to one across the interface within a few grid points. Specifically, the level-set function is converted to the phase-identification function by the following formulas: The above function was chosen because it is differentiable within the entire domain, while it is exactly zero and unity in the fluid and solid domain, respectively. In the volume penalization method, a no-slip condition at the solid surface is realized by introducing an artificial damping force to the Navier-Stokes equation (Equation (5)) as follows: Here, the final term on the right-hand side is the volume penalization term. Obviously, this term has a non-zero value and acts to suppress all the velocity components inside the solid, while Equation (13) reduces to the original Navier-Stokes equation (Equation (5)) within the flow domain where φ = 0. The advantage of the volume penalization method is that solid objects with different shapes can be easily implemented by changing the spatial distribution of φ in the same Cartesian grid system. The drawback is that the grid convergence is relatively slow, since the interface is not explicitly captured and smeared within a few grid points as mentioned above. In the present study, we made grid convergence tests and confirmed that the present results do not change significantly by refining the mesh further ( Figure S1).

Numerical Methods and Conditions
We solved Equations (6) and (13) to obtain the velocity field in the moving coordinate by a pseud-spectral method, in which a solution was expanded by Fourier modes in the x, z direction and Chebyshev polynomials in the y direction, respectively. For time advancement, the Crank-Nicolson was used for diffusion terms, whereas the second-order Adams Bashforth scheme is applied for the convection terms. The inertia and VPM terms appearing in the third and fourth terms are taken into account with the Euler explicit method. The preset code was validated and successfully applied to control and estimation of unsteady turbulent flows in previous studies [24,27].
In this study, we set the diameter and height of micro-pillar to be 200 and 100 µm, respectively. The diameter of the pillar was shown to have only minor effect on the radial profile of the induced flow using the perturbation theory [8], so we studied the flow around the pillar with this representative diameter. The width of the computational domain was 4δ * = 800 µm, which is equal to the spacing of pillars in the experiment. This dimension was set to be large enough so that induced flow profiles of neighboring pillars do not interact. The height of the domain was 2δ * = 200 µm. The numbers of modes employed in the current simulation were (N 1 , N 2 , N 3 ) = (64 × 33 × 64) in x, y, z directions, respectively. 3/2 rule was used for removing aliasing errors, so that the non-linear terms were evaluated in 1.5 times finer physical grid points in each direction. Throughout this work, the vibration amplitude was A = 4 µm and frequency was f = 1000 Hz. Accordingly, dimensionless numbers were Re = 10 and St = 12.4/π, respectively. The numerical time step was set to be ∆t = 1.0 × 10 −4 , which indicates that t −1 = 10 4 time steps are required to compute the velocity field for one oscillation period (t = 1). The computation was started from a stationary flow at t = 0, and the rotational vibration was applied for t = 60 to achieve a fully developed velocity field. After the transient period, the flow field became completely periodic in one oscillation cycle. All statistics shown below were obtained by integrating the velocity data over one oscillation period after the flow field had reached the statistically steady state.

Derivation of Steady Streaming Flow (Time-Averaged Velocity Field)
It is widely known that the SS flow field could be essentially different depending on whether the averaging is made at a fixed stationary location (Eulerian frame) or along a particle moving with the local fluid velocity (Lagrangian frame) due to its oscillatory nature [20,23]. Therefore, we examined and compared the time-averaged SS velocity fields obtained in the Eulerian and Langrangian frames. In the Eulerian approach, the average velocity field was calculated by simply averaging the vector at identical positions in the stationary coordinate system converted from the moving coordinate system.
In the Lagrangian approach, trajectories of virtual fluid particles initially located at a uniform spacing of ∆ = 8.33 µm in the 3D flow field were calculated using the 4th order Runge-Kutta method [28]. The instantaneous velocities between the grids were linearly interpolated from the velocity of the surrounding 6 grid points. The instantaneous velocity field data during one circulating period consisted of 50 time frames, and the velocities at the time points between frames were also linearly interpolated using two neighboring frames. After tracking for five periods of the rotational vibration, the velocity field was obtained from the displacement vectors that connect the start and end points at the identical phase.

Fabrication of Micro-Pillar Array
We fabricate the 5 × 5 array of cylindrical pillars with 200 µm diameter and 100 µm height, with 800 µm center-to-center intervals in accordance with the numerical simulation described in Section 2, using poly-dimethylsiloxane (PDMS) as the material ( Figure 2). Four cylindrical pillars with 1.4 mm diameter and 200 µm height were arranged at the four corners of the substrate as spacers to determine the height of the fluid volume. In practice, the master mold was fabricated on a 2-inch silicon wafer by the deep reactive ion etching apparatus (RIE-400iPB, Samco, Japan) using the Bosch process at a rate of 0.4 µm/cycle. PDMS resin (KE-106, Shin-Etsu Chemical, Japan) mixed with its curing agent at 10:1 weight ratio was poured into the master mold. After curing at 50 • C for 120 min, the PDMS substrate was obtained by peeling it off from the mold (Figure 2b).
It is widely known that the SS flow field could be essentially different depending on whether the averaging is made at a fixed stationary location (Eulerian frame) or along a particle moving with the local fluid velocity (Lagrangian frame) due to its oscillatory nature [20,23]. Therefore, we examined and compared the time-averaged SS velocity fields obtained in the Eulerian and Langrangian frames. In the Eulerian approach, the average velocity field was calculated by simply averaging the vector at identical positions in the stationary coordinate system converted from the moving coordinate system. In the Lagrangian approach, trajectories of virtual fluid particles initially located at a uniform spacing of Δ = 8.33 µm in the 3D flow field were calculated using the 4th order Runge-Kutta method [28]. The instantaneous velocities between the grids were linearly interpolated from the velocity of the surrounding 6 grid points. The instantaneous velocity field data during one circulating period consisted of 50 time frames, and the velocities at the time points between frames were also linearly interpolated using two neighboring frames. After tracking for five periods of the rotational vibration, the velocity field was obtained from the displacement vectors that connect the start and end points at the identical phase.

Fabrication of Micro-Pillar Array
We fabricate the 5 × 5 array of cylindrical pillars with 200 µm diameter and 100 µm height, with 800 µm center-to-center intervals in accordance with the numerical simulation described in Section 2, using poly-dimethylsiloxane (PDMS) as the material (Figure 2). Four cylindrical pillars with 1.4 mm diameter and 200 µm height were arranged at the four corners of the substrate as spacers to determine the height of the fluid volume. In practice, the master mold was fabricated on a 2-inch silicon wafer by the deep reactive ion etching apparatus (RIE-400iPB, Samco, Japan) using the Bosch process at a rate of 0.4 µm/cycle. PDMS resin (KE-106, Shin-Etsu Chemical, Japan) mixed with its curing agent at 10:1 weight ratio was poured into the master mold. After curing at 50 °C for 120 min, the PDMS substrate was obtained by peeling it off from the mold (Figure 2b).

Experimental Setup and Conditions
Since the PDMS is hydrophobic, air bubbles are often trapped around pillars when the liquid (water) is dropped on its surface. Thus, the substrate was made hydrophilic by the oxygen plasma (SEDE-GE, Meiwafosis Co., Ltd., Japan) treatment for 5 min. Immediately after this treatment, 10 µL of deionized (DI) water containing 0.5 µm yellow-green fluorescent beads (F8813, Thermo Fisher Scientific Inc., MA, USA) as a tracer was dropped to the center of the substrate. Then, the substrate was covered with a cover glass. The thickness of the fluid layer was set to 200 µm by spacers (Figure 3a). Then this assembly was fixed to the XY piezo-drive stage (ML-20XYL, MESS-TEK, Japan) ( Figure 3b) using small slips of double-sided tape.
piezo stage using the waveform generator (AG 1022F, OWON, China) via the amplifier (M 2501-1, MESS-TEK, Japan). Applied voltage at 60 V induced A = 4 µm displacement over a wide range of frequency below the resonance of this actuator, which was confirmed from the image of the highspeed camera (Mi-2000, Photron, Japan). The vibration frequency was set at f = 1000 Hz throughout the present study in accordance with the numerical simulation.

PIV System
We used the confocal micro-PIV technique to measure the flow field around a micro-pillar. Figure 4 shows the schematic diagram of the confocal micro-PIV system. In this system, sequential images of fluorescence tracer particles are obtained by a high-speed camera (Mi-2000, Photron, Japan) via the high NA objective water immersion lens (XLUMPLFLN 20 XW, OLYMPUS, Japan) and a confocal scanner (CSU-X1, YOKOGAWA, Japan), and are stored in the PC. The continuous wave (CW) blue laser (488 nm, Sapphire SF, COHERENT, CA, USA) was used as the illumination light source. Since the frequency of the micro-pillar was 1000 Hz, the shutter speed was set to 1/21000 s (~48 µs) so that the instantaneous particle images could be resolved without blurring. The frame rate of the image acquisition was set to 2000 fps. To generate the circular vibration, sinusoidal wave signals with 90 • phase offset was applied to the piezo stage using the waveform generator (AG 1022F, OWON, China) via the amplifier (M 2501-1, MESS-TEK, Japan). Applied voltage at 60 V induced A = 4 µm displacement over a wide range of frequency below the resonance of this actuator, which was confirmed from the image of the high-speed camera (Mi-2000, Photron, Japan). The vibration frequency was set at f = 1000 Hz throughout the present study in accordance with the numerical simulation.

PIV System
We used the confocal micro-PIV technique to measure the flow field around a micro-pillar. Figure 4 shows the schematic diagram of the confocal micro-PIV system. In this system, sequential images of fluorescence tracer particles are obtained by a high-speed camera (Mi-2000, Photron, Japan) via the high NA objective water immersion lens (XLUMPLFLN 20 XW, OLYMPUS, Japan) and a confocal scanner (CSU-X1, YOKOGAWA, Japan), and are stored in the PC. The continuous wave (CW) blue laser (488 nm, Sapphire SF, COHERENT, CA, USA) was used as the illumination light source. Since the frequency of the micro-pillar was 1000 Hz, the shutter speed was set to 1/21,000 s (~48 µs) so that the instantaneous particle images could be resolved without blurring. The frame rate of the image acquisition was set to 2000 fps.

PIV Analysis
Based on the acquired images, we obtained a two-dimensional velocity field using PIV analysis software (Koncerto II, Seika Digital Image, Japan). For obtaining the velocity vector, a recursive crosscorrelation method was used, with a 8 × 8 pixel interrogation window and 50% overlap. This window

PIV Analysis
Based on the acquired images, we obtained a two-dimensional velocity field using PIV analysis software (Koncerto II, Seika Digital Image, Japan). For obtaining the velocity vector, a recursive cross-correlation method was used, with a 8 × 8 pixel interrogation window and 50% overlap. This window size corresponds to 6.4 µm × 6.4 µm in the physical dimension. In the post processing, standard deviation validation and median filter were used to remove incorrect vectors. Since the micro-pillar oscillates at 1000 Hz, two images were recorded during the one cycle of a pillar rotation when the frame rate was 2000 fps. The SS velocity field was obtained from the displacement of tracer images at every two frames; i.e., images at identical rotational phase. Finally, the average velocity field induced around the pillar was obtained by averaging the flow fields of the 60 rotational cycles. We confirmed that the resultant average profile was almost independent of the interrogation window size ( Figure S2).

Horizontal Visualization
A horizontal view of the trajectory of tracer particle was obtained through the objective lens with a long working distance (PAL-10-A, x10, SIGMAKOKI CO. LTD., Japan; WD = 34 mm) located on the side of the pillar substrate and captured by a monochrome CCD camera (BFS-U3-32S4M-C, FLIR Systems, Inc., OR, USA). Polystyrene beads with 10 µm diameter (01-00-104, micromod Partikeltechnologie GmbH., Germany) were used as tracer particles for easy visualization with low magnification lens. The sedimentation velocity of 10 µm bead with a specific density of 1.03 in water is estimated to be around 6.5 µm/s, so that the sedimentation distance within 1.7 s (duration of observation) is not significant.

Instantaneous Velocity Field
In the present numerical simulation, a 3D instantaneous flow field around the pillar was obtained after converting the results in the moving coordinate to the stationary one. Figure 5 shows a two-dimensional (2D) vector plot in the plane 0.5δ away from the non-slip boundary at the bottom (y = 50 µm; the center plane of the pillar in y direction). The vectors and color maps in Figure 5 represent the flow direction and the absolute values of the horizontal component of the velocity (|V hor | = u 2 1 + u 2 3 ), respectively. The region near the pillar is enlarged while the calculation area is 800 µm for x and z directions. Figure 5a-d show the instantaneous fields with the phase shift of every 90 • (see Movie S1 for all velocity fields within one rotation). Each flow field is rotationally symmetric, showing that the flow field reaches a fully developed state. High-speed regions appeared at the front and back faces of the cylindrical pillar in the traveling direction. These high-speed regions separate near the two sides of the pillar and cause a pair of vortex-like structures. The maximum fluid velocity in this region is~25 mm/s, which is comparable to the velocity of oscillation U max = 2πAf 25.2 mm/s. Although not shown in the figure, the comparable speed was also observed at the vicinity of no-slip boundaries (e.g., the upper and lower walls). and back faces of the cylindrical pillar in the traveling direction. These high-speed regions separate near the two sides of the pillar and cause a pair of vortex-like structures. The maximum fluid velocity in this region is ~ 25 mm/s, which is comparable to the velocity of oscillation Umax = 2πAf ~ 25.2 mm/s. Although not shown in the figure, the comparable speed was also observed at the vicinity of no-slip boundaries (e.g., the upper and lower walls).

Time-Averaged Flow Field
In a periodically vibrating flow field, the time-averaged velocity becomes null if there is no obstacle in the oscillating direction. However, a non-zero net velocity field appears when the obstacle exists as a result of the Stokes drift. As mentioned previously, we examined the two time-averaging approaches, i.e., Eulerian and Lagrangian averaging. In the latter, we calculated the velocity vectors from the displacement of the virtual tracer particles. An example of the 2D trajectory of a tracer, initially placed at (x, y, z) = (150, 50, 0 µm) (horizontally 50 µm away from the side wall of the pillar) is shown in Figure 6. The tracer particle moves along the distorted orbital path, which can be seen as the superposition of the circular periodic motion and the steady translational movement toward the upper left in the figure. When the pillar is absent, the trajectory draws a perfect circle and returns to the original position after one rotation. Due to the existence of the pillar, however, the position in the same phase is shifted after each rotation. This displacement divided by the time of one rotation period corresponds to the mean velocity observed in experimental particle tracking generated in the SS flow field.

Time-Averaged Flow Field
In a periodically vibrating flow field, the time-averaged velocity becomes null if there is no obstacle in the oscillating direction. However, a non-zero net velocity field appears when the obstacle exists as a result of the Stokes drift. As mentioned previously, we examined the two time-averaging approaches, i.e., Eulerian and Lagrangian averaging. In the latter, we calculated the velocity vectors from the displacement of the virtual tracer particles. An example of the 2D trajectory of a tracer, initially placed at (x, y, z) = (150, 50, 0 µm) (horizontally 50 µm away from the side wall of the pillar) is shown in Figure 6. The tracer particle moves along the distorted orbital path, which can be seen as the superposition of the circular periodic motion and the steady translational movement toward the upper left in the figure. When the pillar is absent, the trajectory draws a perfect circle and returns to the original position after one rotation. Due to the existence of the pillar, however, the position in the same phase is shifted after each rotation. This displacement divided by the time of one rotation period corresponds to the mean velocity observed in experimental particle tracking generated in the SS flow field.
The two-dimensional vector plots of the averaged velocity field obtained by Eulerian and Lagrangian methods at the y = 50 µm horizontal plane are depicted in Figure 7. The color map represents the absolute value of the horizontal component |V hor | of the average velocity. It is clear that a net flow is induced around the pillar in both cases. However, the radial peak position of the velocity is closer to the pillar in the former case (r~120 µm in the Eulerian method, but r~130 µm in the Lagrangian method). Moreover, the peak value is about two times greater in Eulerian method compared to the Lagrangian method. The result clearly indicates that the averaged flow field depends on the averaging methods. Because the mass transport and the paths of suspended molecules/particles are governed by the Lagrangian trajectories, the Lagrangian averaging is necessary to predict the above-mentioned transport phenomena in micro devices. In the following sections, the Lagrangian velocity obtained from the present simulation will be compared with the averaged translational velocity of fluorescence particles in the experiment in order to validate the present numerical code. The two-dimensional vector plots of the averaged velocity field obtained by Eulerian and Lagrangian methods at the y = 50 µm horizontal plane are depicted in Figure 7. The color map represents the absolute value of the horizontal component |Vhor| of the average velocity. It is clear that a net flow is induced around the pillar in both cases. However, the radial peak position of the velocity is closer to the pillar in the former case (r ~ 120 µm in the Eulerian method, but r ~ 130 µm in the Lagrangian method). Moreover, the peak value is about two times greater in Eulerian method compared to the Lagrangian method. The result clearly indicates that the averaged flow field depends on the averaging methods. Because the mass transport and the paths of suspended molecules/particles are governed by the Lagrangian trajectories, the Lagrangian averaging is necessary to predict the above-mentioned transport phenomena in micro devices. In the following sections, the Lagrangian velocity obtained from the present simulation will be compared with the averaged translational velocity of fluorescence particles in the experiment in order to validate the present numerical code.   The two-dimensional vector plots of the averaged velocity field obtained by Eulerian and Lagrangian methods at the y = 50 µm horizontal plane are depicted in Figure 7. The color map represents the absolute value of the horizontal component |Vhor| of the average velocity. It is clear that a net flow is induced around the pillar in both cases. However, the radial peak position of the velocity is closer to the pillar in the former case (r ~ 120 µm in the Eulerian method, but r ~ 130 µm in the Lagrangian method). Moreover, the peak value is about two times greater in Eulerian method compared to the Lagrangian method. The result clearly indicates that the averaged flow field depends on the averaging methods. Because the mass transport and the paths of suspended molecules/particles are governed by the Lagrangian trajectories, the Lagrangian averaging is necessary to predict the above-mentioned transport phenomena in micro devices. In the following sections, the Lagrangian velocity obtained from the present simulation will be compared with the averaged translational velocity of fluorescence particles in the experiment in order to validate the present numerical code.

Results of PIV Measurement
The image of the movement of fluorescent tracer particles obtained at the y = 50 µm horizontal plane is shown in Figure 8a. Sixty instantaneous images at an identical rotational phase, captured by a high-speed camera, are superimposed. Although tracer beads actually moved along the rotating trajectory in accordance with the rotational vibration, the net displacement along the pillar sidewall can be clearly seen by connecting the positions at the same phase. The result of PIV analysis, obtained in this series of images, is depicted in Figure 8b. The condition is the same as that in the numerical result shown in Figure 7. The overall profile and the magnitude of the velocity are similar to those obtained in the Lagrangian averaging (Figure 7b). trajectory in accordance with the rotational vibration, the net displacement along the pillar sidewall can be clearly seen by connecting the positions at the same phase. The result of PIV analysis, obtained in this series of images, is depicted in Figure 8b. The condition is the same as that in the numerical result shown in Figure 7. The overall profile and the magnitude of the velocity are similar to those obtained in the Lagrangian averaging (Figure 7b).

Comparison of Radial Velocity Profile
To compare the results quantitatively, we plotted the radial profiles of |Vhor| in Figure 9. The profiles were averaged in the azimuthal direction of the pillar. The magnitude of the peak velocity obtained by Eulerian method is twice as large as the other, and the peak position of the velocity is closer to the pillar, as qualitatively seen in Figure 7. On the other hand, the velocity distributions obtained by the Lagrangian method and PIV measurement showed a similar trend in terms of the peak position (r = 130 µm; 30 µm away from the pillar wall) and the decay of the profile. There is a difference between the two profiles close to the pillar; this difference could be partly caused by the difficulty in resolving the high-shear region in both simulation and PIV. The size of the single grid in the simulation corresponds to be 8.3 µm, and the size of the window for image correlation in PIV corresponds to be 6.4 µm. Particle images tended to blur in the region close to the wall due to the 3D flow described later. Nonetheless, a good agreement in the decaying profile after the peak supports the validity of the present simulation.
Overall, the significant difference between two averaging methods is reasonable considering that the flow is unsteady. The virtual fluid particles in the flow do not stay at the constant radial position r with respect to the center of the pillar. Instead, their radial position changes as they rotated around the reference position (center of the rotational path) at each phase. Thus, the velocity at the actual particle position is different from that at the reference position. As this slight difference accumulates, the time averaged velocity per one cycle becomes different from the average velocity obtained at the stationary coordinates. The present result confirms that it is necessary to use the Lagrangian tracking method to reproduce the unsteady flow around the vibrating pillar.

Comparison of Radial Velocity Profile
To compare the results quantitatively, we plotted the radial profiles of |V hor | in Figure 9. The profiles were averaged in the azimuthal direction of the pillar. The magnitude of the peak velocity obtained by Eulerian method is twice as large as the other, and the peak position of the velocity is closer to the pillar, as qualitatively seen in Figure 7. On the other hand, the velocity distributions obtained by the Lagrangian method and PIV measurement showed a similar trend in terms of the peak position (r = 130 µm; 30 µm away from the pillar wall) and the decay of the profile. There is a difference between the two profiles close to the pillar; this difference could be partly caused by the difficulty in resolving the high-shear region in both simulation and PIV. The size of the single grid in the simulation corresponds to be 8.3 µm, and the size of the window for image correlation in PIV corresponds to be 6.4 µm. Particle images tended to blur in the region close to the wall due to the 3D flow described later. Nonetheless, a good agreement in the decaying profile after the peak supports the validity of the present simulation.

Comparison of Normal Velocity Profile
Next, we examined the velocity distribution in the vertical plane including the Stokes layer. To do this we conducted the micro PIV measurement from y = 0 to 100 µm at 10 µm interval. The plane

Comparison of Normal Velocity Profile
Next, we examined the velocity distribution in the vertical plane including the Stokes layer. To do this we conducted the micro PIV measurement from y = 0 to 100 μm at 10 μm interval. The plane ) PIV measurement. PIV was repeated three times with different setup, and the average and standard deviation are shown.
Overall, the significant difference between two averaging methods is reasonable considering that the flow is unsteady. The virtual fluid particles in the flow do not stay at the constant radial position r with respect to the center of the pillar. Instead, their radial position changes as they rotated around the reference position (center of the rotational path) at each phase. Thus, the velocity at the actual particle position is different from that at the reference position. As this slight difference accumulates, the time averaged velocity per one cycle becomes different from the average velocity obtained at the stationary coordinates. The present result confirms that it is necessary to use the Lagrangian tracking method to reproduce the unsteady flow around the vibrating pillar.

Comparison of Normal Velocity Profile
Next, we examined the velocity distribution in the vertical plane including the Stokes layer. To do this we conducted the micro PIV measurement from y = 0 to 100 µm at 10 µm interval. The plane at y > 100 µm was not obtained experimentally due to the limitation in the working distance of the objective. The color maps of |V hor | in r-y plane obtained in the numerical simulation (Lagrangian averaging) and the PIV measurement are shown in Figure 10a,b, respectively. They showed a similar trend in terms of the peak position of r = 130 µm at 40 < y < 80 µm, which slightly moves closer to the pillar near the top (y~100 µm). This is attributed to the strong three-dimensionality of the flow around the tip of the pillar, as will be shown later. The magnitude of the velocity decays rapidly as the position gets closer to the bottom wall.

Comparison of Normal Velocity Profile
Next, we examined the velocity distribution in the vertical plane including the Stokes layer. To do this we conducted the micro PIV measurement from y = 0 to 100 µm at 10 µm interval. The plane at y > 100 µm was not obtained experimentally due to the limitation in the working distance of the objective. The color maps of |Vhor| in r-y plane obtained in the numerical simulation (Lagrangian averaging) and the PIV measurement are shown in Figure 10a,b, respectively. They showed a similar trend in terms of the peak position of r = 130 µm at 40 < y < 80 µm, which slightly moves closer to the pillar near the top (y ~ 100 µm). This is attributed to the strong three-dimensionality of the flow around the tip of the pillar, as will be shown later. The magnitude of the velocity decays rapidly as the position gets closer to the bottom wall. The horizontal velocity distributions in the vertical (y-) direction at r = 130 µm are plotted shown in Figure 11. Both profiles agreed well, with velocity values close to 0 µm/s on the substrate surface The horizontal velocity distributions in the vertical (y-) direction at r = 130 µm are plotted shown in Figure 11. Both profiles agreed well, with velocity values close to 0 µm/s on the substrate surface (y = 0 µm), which increased to~1.8 mm/s towards the tip of the pillar. The peak value is obtained in the vicinity of y = 70 to 80 µm.
In the Stokes second problem of a unidirectional oscillation of the flat solid wall in contact with the fluid, the fluid in the vicinity of the wall almost follows the wall movement, but this influence exponentially decays as it moves away from the wall with a characteristic length δ s = (2ν/ω) 1/2 [13]. From its analytical solution, the effect of the wall motion decreases down to~5% (1/e 3 ) at y = 3δ s . In the present condition, 3δ s corresponds to 53.4 µm by substituting ν = 1.0 × 10 −6 m 2 /s (25 • C) and ω = 2πf = 6.28 × 10 3 rad/s. The height is indicated by dotted lines in Figures 10 and 11. Since the net flow discussed in the present system is induced by the relative velocity between the fluid and the substrate, the maximum average flow field occurs outside the Stokes layer. In contrast, the fluid within the Stokes layer thickness 3δ s is dragged by the motion of the wall due to the viscous effect so that the relative velocity diminishes on the bottom wall.
The present results indicate that the viscous effect adjacent to the solid wall has a strong impact on the net travel distance of a passive tracer, and the thickness of this fluid layer is determined by the Stokes layer thickness. It should be also noted that, for each height y, the peak of the averaged velocity occurs about 30 µm away from the side wall of the pillar as shown in Figure 10. This horizontal gap is also similar to the Stokes layer thickness. When SS is applied to micro devices, the Stokes layer thickness is not negligibly small in general, and this causes the strong dependency of the induced averaged flow on the vertical location, i.e., the distance from the top and bottom substrate. Therefore, in order to predict and design the entire velocity field inside a micro device, 3D analysis is requisite.
substrate, the maximum average flow field occurs outside the Stokes layer. In contrast, the fluid within the Stokes layer thickness 3δs is dragged by the motion of the wall due to the viscous effect so that the relative velocity diminishes on the bottom wall.
The present results indicate that the viscous effect adjacent to the solid wall has a strong impact on the net travel distance of a passive tracer, and the thickness of this fluid layer is determined by the Stokes layer thickness. It should be also noted that, for each height y, the peak of the averaged velocity occurs about 30 µm away from the side wall of the pillar as shown in Figure 10. This horizontal gap is also similar to the Stokes layer thickness. When SS is applied to micro devices, the Stokes layer thickness is not negligibly small in general, and this causes the strong dependency of the induced averaged flow on the vertical location, i.e., the distance from the top and bottom substrate. Therefore, in order to predict and design the entire velocity field inside a micro device, 3D analysis is requisite. Figure 11. Comparison of vertical distributions of mean horizontal velocity magnitude |Vhor|.

Three-Dimensionality of the Flow
Lastly, we investigated the three-dimensionality of the averaged velocity field around the pillar. Figure 12 shows a vector and a contour plot of the averaged velocity in the vertical plane (| | = + ) obtained from the numerical simulation. In the vicinity of the upper and lower wall surfaces, it is clear that there is almost no upward/downward motion, so that the flow is practically two-dimensional. However, at the intermediate region between the top and bottom plates, the velocity toward the pillar with the magnitude as large as ~200 µm/s appears around the corner of the pillar (r = 130 to 150 µm and y = 50 to 100 µm). This flow is diverted upward at the pillar edge. As Figure 11. Comparison of vertical distributions of mean horizontal velocity magnitude |V hor |.

Three-Dimensionality of the Flow
Lastly, we investigated the three-dimensionality of the averaged velocity field around the pillar. Figure 12 shows a vector and a contour plot of the averaged velocity in the vertical plane (|V ver | = u 2 1 + u 2 2 ) obtained from the numerical simulation. In the vicinity of the upper and lower wall surfaces, it is clear that there is almost no upward/downward motion, so that the flow is practically two-dimensional. However, at the intermediate region between the top and bottom plates, the velocity toward the pillar with the magnitude as large as~200 µm/s appears around the corner of the pillar (r = 130 to 150 µm and y = 50 to 100 µm). This flow is diverted upward at the pillar edge. As a result, a vortical motion is generated close to the corner of the pillar. Although the magnitude is about 10 times smaller than |V hor |, the presence of the vertical velocity component has a strong impact on the averaged velocity field as discussed in Figure 10.
To confirm the three-dimensionality of the averaged velocity flow field in the experimental flow field, the motion of the particle was observed from the side of the pillar using an objective lens with a long working distance placed horizontally. Because the confocal illumination is not possible away from the wall, we introduced 10 µm polystyrene beads as tracers into the fluid at a low concentration so that individual beads were visible from a long distance. Two representative paths visualized by superimposing~50 successive frames (1.7 s) are shown on the left column in Figure 13a,b (see Movie S2 and S3 for corresponding movies). In both cases, tracer particles exhibited three-dimensional motions including upward and downward movements, instead of a simple orbital movement in the same horizontal plane. In Figure 13a, a bead initially circulating around the middle of the pillar was suddenly raised at the region close the apex, and then hovered above the pillar. The similar path could be reproduced by tracking an ideal tracer in the averaged velocity field obtained by the Lagrangian averaging (Figures 7b and 12) in the simulation (right figures). This sudden ascension should be caused by the upward flow near the apex discussed in Figure 12. Another bead shown in Figure 13b, which was circulating a bit far outside of the pillar slightly, shifted downward as it got closer to the pillar, and then raised. This path was also reproduced in the numerical simulation (right figures), which was convected by the downward net velocity at r = 120~170 µm and y = 100~130 µm, and raised by the upward net velocity at r = 110~150 µm and y = 40~100 µm shown in Figure 12. These trajectories show the validity of numerical prediction for three-dimensionality of flow induced around the pillar at the present vibration condition.
Micromachines 2018, 9, x FOR PEER REVIEW 15 of 18 a result, a vortical motion is generated close to the corner of the pillar. Although the magnitude is about 10 times smaller than |Vhor|, the presence of the vertical velocity component has a strong impact on the averaged velocity field as discussed in Figure 10. To confirm the three-dimensionality of the averaged velocity flow field in the experimental flow field, the motion of the particle was observed from the side of the pillar using an objective lens with a long working distance placed horizontally. Because the confocal illumination is not possible away from the wall, we introduced 10 µm polystyrene beads as tracers into the fluid at a low concentration so that individual beads were visible from a long distance. Two representative paths visualized by superimposing ~50 successive frames (1.7 s) are shown on the left column in Figure 13a,b (see Movie S2 and S3 for corresponding movies). In both cases, tracer particles exhibited three-dimensional motions including upward and downward movements, instead of a simple orbital movement in the same horizontal plane. In Figure 13a, a bead initially circulating around the middle of the pillar was suddenly raised at the region close the apex, and then hovered above the pillar. The similar path could be reproduced by tracking an ideal tracer in the averaged velocity field obtained by the Lagrangian averaging (Figures 7b and 12) in the simulation (right figures). This sudden ascension should be caused by the upward flow near the apex discussed in Figure 12. Another bead shown in Figure 13b, which was circulating a bit far outside of the pillar slightly, shifted downward as it got closer to the pillar, and then raised. This path was also reproduced in the numerical simulation (right figures), which was convected by the downward net velocity at r = 120 ~ 170 µm and y = 100 ~ 130 µm, and raised by the upward net velocity at r = 110 ~ 150 µm and y = 40 ~ 100 µm shown in Figure 12. These trajectories show the validity of numerical prediction for three-dimensionality of flow induced around the pillar at the present vibration condition.

Conclusions
In this work, we developed a simulation tool to predict the SS flow induced by vibration without assuming 2D flow and small vibration amplitude. The developed numerical code was based on the volume penalization method in the Cartesian coordinate system, so that arbitrary three-dimensional structures can be readily embedded without requiring grid generation for each geometry. As pointed out in the past SS studies, we confirmed that the average velocity at fixed points (Eulerian mean) and

Conclusions
In this work, we developed a simulation tool to predict the SS flow induced by vibration without assuming 2D flow and small vibration amplitude. The developed numerical code was based on the volume penalization method in the Cartesian coordinate system, so that arbitrary three-dimensional structures can be readily embedded without requiring grid generation for each geometry. As pointed out in the past SS studies, we confirmed that the average velocity at fixed points (Eulerian mean) and that of fluid particles advected by the local flow velocity (Lagrangian mean) are essentially different. In practical microdevice applications, the Lagrangian-averaged flow governs mass transport and mixing. In this work, we compared the Lagrangian mean field obtained by numerical calculation with micro PIV experimental results. The quantitative agreement between them validates the present simulation method.
Through these numerical and experimental analyses, we observed two 3D characteristics of the SS flow even at the flow Reynolds number of 10. First, inside the Stokes layer, which develops from the solid surface toward the fluid region, the averaged velocity due to the Stokes drift diminishes because the fluid is dragged by the solid motion due to the viscous effect, and thereby the relative velocity between the fluid and the solid reduces. Since the thickness of the Stokes layer cannot be neglected with respect to the length scale of microdevices in general, it is necessary to accurately predict the influence of the Stokes layer in order to understand the flow inside the device. Secondly, a large strong vortex motion accompanying upward and downward motions was confirmed in the vicinity of the apex of the micro-pillar. Such vertical fluid motions could have significant influences on mass transport, and may result in a complicated flow field combined with the three-dimensionality caused by the Stokes layer mentioned above. Since the present numerical results reproduce the 3D nature of the flow fields observed in experiment, our numerical code could be a powerful tool to optimize the structure of a micro-pillar and the vibration mode in order to develop innovative microfluidic devices in future work.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-666X/9/12/668/ s1, Supplementary Text: Governing equations of fluid in a moving frame with rotational vibration. Figure S1: Comparison with the radial profile of |V hor | calculated with different grid resolutions in the numerical simulation. The distribution calculated with a low-resolution grid (48 × 25 × 48) differ significantly, but those calculated with medium (96 × 49 × 96; used in the main results) and high-resolutions (144 × 73 × 144) were similar., Figure S2: Comparison of the radial profile of |V hor | calculated with different window size in PIV analysis. One pixel in images corresponds to 0.8 µm. Supporting Movie 1: Animation of the velocity field (u 1 , u 3 ) at y = 50 µm obtained from the numerical simulation., Supporting Movie 2: Motion of a tracer particle depicted in Figure 13a in the main text., Supporting Movie 3: Motion of a tracer particle depicted in Figure 13b in the main text.