Numerical analysis of flow anisotropy in rotated-square deterministic lateral displacement devices at moderate Reynolds number

Deterministic lateral displacement (DLD) is a microfluidic method for accurately separating particles by size or deformability. Recent efforts to operate DLD devices in the inertial, rather than in the Stokes flow regime, have been hindered by a loss of separation efficiency and difficulty predicting the separation behavior. One factor contributing to these problems is the onset of inertia-induced flow anisotropy where the average flow direction does not align with the direction of the pressure gradient in the device. We use the lattice-Boltzmann method to simulate two-dimensional flow through a rotated-square DLD geometry with circular pillars at a Reynolds number up to 100 for different gap sizes and rotation angles. We find that anisotropy in this geometry is a nonmonotonous function of Reynolds number and can be positive or negative. This finding is in contradiction to the naive expectation that inertia would always drive flow along the principal direction of the pillar array. Anisotropy tends to increase in magnitude with gap size and rotation angle. By analyzing the traction distribution along the pillar surface, we explain how the change of the flow field upon increasing inertia leads to the observed trends of anisotropy. Our work contributes to a better understanding of the inertial flow behavior in ordered cylindrical porous media, and it might contribute to improved DLD designs for operation in the inertial regime.

Figure 1

Illustration of the rotated-square array. The circular pillars of radius $r$ lie on the nodes of a square lattice with center-to-center distance $\lambda$. The gap size between pillars is $G=\lambda -2r$. The array is rotated at an angle $\alpha$ relative to the coordinate axes. By varying the direction of the body force density driving the flow, which is denoted ${\mathit{F}}_{\text{b}}$, the average velocity of the fluid $\overline{\mathit{u}}$ is imposed to be parallel to the basis vector $\widehat{\mathit{x}}$, which indicates the axis of the DLD device. For the validation test presented in Sec. 3a, $\alpha$ is set to zero and the body force ${\mathit{F}}_{\text{b}}$ is maintained at an angle $\beta$. Simulations are carried out over one unit cell of the array with periodic boundary conditions, highlighted by the gray box.

Figure 2

(a) Study of convergence tolerance and grid independence with $G^{★}=0.5, \text{Re}=50$. Doubling the resolution to $400\times 400$ grid points increases the anisotropy by $0.65\%$ at a convergence tolerance of ${10}^{-8}$ between times $t$ and $t-1000\mathrm{\Delta }t$. (b) Average flow direction relative to the array orientation at $G^{★}=0.51$ against $\text{Re}$, for simulations with the flow allowed to freely develop (body force at an angle $\beta$ with respect to $\widehat{\mathit{x}}$). In black are the simulation results obtained by Koch and Ladd [22].

Figure 3

Streamlines colored by the normalized velocity magnitude (see the right color bar) for different $G^{★}$ at $\alpha =5^{\circ }$. At low Re (a)–(c), the flow field is symmetrical and attached to pillar surfaces. At higher Re (d)–(f), a recirculating wake is visible downstream from the pillar. The red line illustrates the path over which the velocity profiles in Fig. 4 are sampled. Note that four periodic unit cells (and therefore four pillars) are shown for visual clarity.

Figure 4

Velocity profiles along the vertical center-to-center line (shown as a red line in Fig. 3) for $\alpha =5^{\circ }$ and different values of $G^{★}$ and Re. The $x$-axis is normalized such that the top side of the pillar corresponds to $-1$ and the bottom side of the neighboring pillar to $+1$. The $y$-axis is normalized by the value of the maximum velocity magnitude.

Figure 5

Anisotropy $A$ for different values of $G^{★}, \alpha$, and Re. Results are grouped by the value of $\alpha$ to highlight the similarity of the results for different values of $\alpha$. Note that the scale of the $y$-axis changes between panels. Points for which the anisotropy is negative are marked with filled shapes.

Figure 6

(a)–(c) Normal and (d)–(f) shear contributions to anisotropy for the simulations shown in Fig. 5. Points for which the anisotropy is below zero are marked with filled shapes.

Figure 7

Polar representation of the normal contribution to the traction (black) and lift (blue) along the pillar surface for (a)–(c) $G^{★}=0.2$, (d)–(f) $G^{★}=0.4$, and (g)–(i) $G^{★}=0.6$ at $\alpha =5^{\circ }$ for different values of Re. The average flow direction is along the $x$-axis. The traction and lift distributions are normalized by their maximum values. Positive values are shown as a solid line, negative values as a dashed line. The value of the normal contribution to the anisotropy is reported in each panel. Panel (a) also contains labels of the quadrants, and the labels for mechanisms associated with each region of the lift distribution, as referred to in Sec. 3d. The normalized pressure field is shown as a color gradient.

Figure 8

Polar representation of the shear contribution to the traction (black) and lift (blue) along the pillar surface for (a)–(c) $G^{★}=0.2$, (d)–(f) $G^{★}=0.4$, and (g)–(i) $G^{★}=0.6$ at $\alpha =5^{\circ }$ for different values of Re. The average flow direction is along the $x$-axis. The traction and lift distributions are normalized by their maximum values. Positive values are shown as a solid line, negative values as a dashed line. The value of the shear contribution to the anisotropy is reported in each panel. The fluid flow is shown as streamlines using the same coloring as in Fig. 3.

Figure 9

Magnitude of normal lift force associated with the different mechanisms identified in Sec. 3d. The sum of mechanisms (i) and (ii) ($m_1$ and $m_2$), associated with the upstream stagnation point and vertical flow and contributing to positive anisotropy, is compared with mechanism (iii) ($m_3$), associated with the skewed horizontal velocity profile and contributing to negative anisotropy. When the solid line is above the broken line, we expect to see an overall negative value of the normal anisotropy, $A^{\text{n}}$. For $G^{★}=0.6$, values are only plotted up to $\text{Re}=80$, as the lift profile has only three zeros at $\text{Re}=100$.

Figure 10

Illustration of the employed interpolation-extrapolation scheme. Lattice points are shown in red (solid nodes) and blue (fluid nodes). Flow information at points 1 and 2 is bilinearly interpolated from the four nodes defining the surrounding voxels, respectively. The interpolated values are then linearly extrapolated to point 0 on the cylinder surface.