Rotation of a low-Reynolds-number watermill: theory and simulations

TL;DR

This study combines theory and simulations to analyze the hydrodynamics of a flow-driven micro-scale watermill, revealing how its rotation depends on the number of rods, with hydrodynamic interactions playing a crucial role.

Contribution

It introduces a combined theoretical and numerical approach to understand the hydrodynamics of a low-Reynolds-number watermill, highlighting the importance of hydrodynamic interactions for accurate predictions.

Findings

Rotational velocity is largely independent of the number of rods for N ≥ 4.

Hydrodynamic interactions significantly affect the watermill's rotation and are essential for accurate modeling.

Theoretical predictions align with simulations when hydrodynamic interactions are included.

Abstract

Recent experiments have demonstrated that small-scale rotary devices installed in a microfluidic channel can be driven passively by the underlying flow alone without resorting to conventionally applied magnetic or electric fields. In this work, we conduct a theoretical and numerical study on such a flow-driven "watermill" at low Reynolds number, focusing on its hydrodynamic features. We model the watermill by a collection of equally-spaced rigid rods. Based on the classical resistive force (RF) theory and direct numerical simulations, we compute the watermill's instantaneous rotational velocity as a function of its rod number $N$, position and orientation. When $N \geq 4$, the RF theory predicts that the watermill's rotational velocity is independent of $N$ and its orientation, implying the full rotational symmetry (of infinity order), even though the geometrical configuration exhibits a lower-fold rotational symmetry; the numerical solutions including hydrodynamic interactions show a weak dependence on $N$ and the orientation. In addition, we adopt a dynamical system approach to identify the equilibrium positions of the watermill and analyse their stability. We further compare the theoretically and numerically derived rotational velocities, which agree with each other in general, while considerable discrepancy arises in certain configurations owing to the hydrodynamic interactions neglected by the RF theory. We confirm this conclusion by employing the RF-based asymptotic framework incorporating hydrodynamic interactions for a simpler watermill consisting of two or three rods and we show that accounting for hydrodynamic interactions can significantly enhance the accuracy of the theoretical predictions.

Tables

Table 1. Table 1: The nondimensional hydrodynamic torque T ¯ z = T z / ( μ ​ ω ​ L 3 ) subscript ¯ 𝑇 𝑧 subscript 𝑇 𝑧 𝜇 𝜔 superscript 𝐿 3 \bar{T}_{z}=T_{z}/\left(\mu\omega L^{3}\right) exerted on a single rod of varying slenderness ϵ italic-ϵ \epsilon which rotates about one of its ends at a prescribed rotational velocity ω 𝜔 \omega in the x ​ y 𝑥 𝑦 xy plane.

pread 0pt —X[3l]?X[c]—X[c]—X[c]—X[c]—X[c]—	$ϵ=0.01$	$ϵ=0.02$	$ϵ=0.03$	$ϵ=0.04$	$ϵ=0.05$
$T_z/\mu \omega L^3$ (RF theory)	0.9096	1.0707	1.1946	1.3013	1.3983
$T_z/\mu \omega L^3$ (Simulation)	0.9093	1.0758	1.1998	1.3088	1.3982

Table 2. Table 2: The nondimensional rotational velocity Ω z / γ ˙ subscript Ω 𝑧 ˙ 𝛾 \Omega_{z}/\dot{\gamma} of a prolate ( κ = 2 𝜅 2 \kappa=2 ) and an oblate ( κ = 0.6 𝜅 0.6 \kappa=0.6 ) spheroid in shear flow 𝐔 ∞ = − γ ˙ ​ y ​ 𝐞 x superscript 𝐔 ˙ 𝛾 𝑦 subscript 𝐞 𝑥 \mathbf{U}^{\infty}=-\dot{\gamma}y\mathbf{e}_{x} as a function of the angle α 𝛼 \alpha between its axis of revolution and the streamwise ( x 𝑥 x ) direction. The theoretical results are from Jeffery ( 1922 ) .

pread 0pt —X[2l]?X[c]X[c]X[c]X[c]— $\kappa =2$ (Prolate)	$\alpha =0^{\circ }$	$\alpha ={30}^{\circ }$	$\alpha ={60}^{\circ }$	$\alpha ={90}^{\circ }$
${\Omega }_z/\dot{\gamma }$ (Theory)	0.1111	0.3056	0.6944	0.8889
${\Omega }_z/\dot{\gamma }$ (Simulation)	0.1055	0.3025	0.6974	0.8945
$\kappa =0.6$ (Oblate)	$\alpha =0^{\circ }$	$\alpha ={30}^{\circ }$	$\alpha ={60}^{\circ }$	$\alpha ={90}^{\circ }$
${\Omega }_z/\dot{\gamma }$ (Theory)	0.8224	0.6612	0.3388	0.1776
${\Omega }_z/\dot{\gamma }$ (Simulation)	0.8293	0.6646	0.3354	0.1705

Equations

$$\beta$$

$$\rho$$

$${\mathbf{f}}^{(k)}(s)=-{\mathbf{J}}^{(k)}\cdot \left[\frac{\partial {\mathbf{r}}^{(k)}}{\partial t}-{\mathbf{U}}^{\infty }\left({\mathbf{r}}^{(k)}\right)\right],$$

$${\bar{T}}_z^{(k)}=\frac{1}{3{\rho }^2\ln {ϵ}^{-1}}\left[6\left({\beta }^2-{\rho }^2\right)\cos {\theta }^{(k)}+8\beta {\cos }^2{\theta }^{(k)}+3{\cos }^3{\theta }^{(k)}+4{\rho }^2{\bar{\omega }}_z\right],$$

$$\displaystyle\displaystyle 3\rho^{2}\left(\ln{{\epsilon}^{-1}}\right)\bar{T}_{z}=\sum_{k=1}^{N}\Bigg{\{}\left[6\left(\beta^{2}-\rho^{2}\right)+\frac{9}{4}\right]\cos{\theta^{(k)}}+$$

$$\displaystyle 4\beta\cos{2\theta^{(k)}}+\frac{3}{4}\cos{3\theta^{(k)}}\Bigg{\}}+4N\left(\beta+\rho^{2}\bar{\omega}_{z}\right).$$

$$\sum _{k=1}^N\cos {\theta }^{(k)}=\{\begin{array}{ll}\cos \theta ,& N=1\\ 0,& N\ge 2,\end{array}$$

$$\sum _{k=1}^N\cos 2{\theta }^{(k)}=\{\begin{array}{ll}\cos 2\theta ,& N=1\\ 2\cos 2\theta ,& N=2\\ 0,& N\ge 3,\end{array}$$

$$\sum _{k=1}^N\cos 3{\theta }^{(k)}=\{\begin{array}{ll}\cos 3\theta ,& N=1\\ 0,& N=2\\ 3\cos 3\theta ,& N=3\\ 0,& N\ge 4,\end{array}$$

$$3{\rho }^2\left(\ln {ϵ}^{-1}\right){\bar{T}}_z=\{\begin{array}{ll}[6({\beta }^2-{\rho }^2)+\frac{9}{4}]\cos \theta +4\beta \cos 2\theta +\frac{3}{4}\cos 3\theta +4(\beta +{\rho }^2{\bar{\omega }}_z),& N=1\\ 8\beta \cos 2\theta +8(\beta +{\rho }^2{\bar{\omega }}_z),& N=2\\ \frac{9}{4}\cos 3\theta +12(\beta +{\rho }^2{\bar{\omega }}_z),& N=3\\ 4N(\beta +{\rho }^2{\bar{\omega }}_z),& N\ge 4.\end{array}$$

$${\bar{\omega }}_z=\{\begin{array}{ll}-\cos \theta [6({\beta }^2-{\rho }^2)+8\beta \cos \theta +3{\cos }^2\theta ]/\left(4{\rho }^2\right),& N=1\\ -2\beta {\cos }^2\theta /{\rho }^2,& N=2\\ -\left[16\beta +3\cos 3\theta \right]/\left(16{\rho }^2\right),& N=3\\ -\beta /{\rho }^2,& N\ge 4.\end{array}$$

$$-\nabla \bar{p}+{\nabla }^2\bar{\mathbf{u}}$$

$$\nabla \cdot \bar{\mathbf{u}}$$

$$\bar{\mathbf{u}}\left({\bar{\mathbf{x}}}_w\right)={\bar{\omega }}_z{\mathbf{e}}_z\times ({\bar{\mathbf{x}}}_w-{\bar{\mathbf{x}}}_{\mathrm{joint}}),$$

$${\int }_{S_w}({\bar{\mathbf{x}}}_w-{\bar{\mathbf{x}}}_{\mathrm{joint}})\times \left[\left(-\bar{p}\mathbf{I}+\nabla \bar{\mathbf{u}}+\nabla {\bar{\mathbf{u}}}^T\right)\cdot {\mathbf{n}}_w\right]\mathrm{d}S=0,$$

$${\mathbf{f}}^{(1)}-{\mathbf{J}}^{(1)}\cdot {\mathbf{v}}^{(2)\to (1)}=-{\mathbf{J}}^{(1)}\cdot \left[\frac{\partial {\mathbf{r}}^{(1)}}{\partial t}-{\mathbf{U}}^{\infty }({\mathbf{r}}^{(1)})\right],$$

$${\mathbf{f}}^{(2)}-{\mathbf{J}}^{(2)}\cdot {\mathbf{v}}^{(1)\to (2)}=-{\mathbf{J}}^{(2)}\cdot \left[\frac{\partial {\mathbf{r}}^{(2)}}{\partial t}-{\mathbf{U}}^{\infty }({\mathbf{r}}^{(2)})\right],$$

$${\mathbf{v}}^{(2)\to (1)}=\frac{1}{4\pi \mu }\ln \left(\frac{h(s)}{L}\right)(\mathbf{I}+{\mathbf{e}}_y{\mathbf{e}}_y)\cdot {\mathbf{f}}^{(2)}(s),$$

$$-{\mathbf{J}}^{(1)}\cdot {\mathbf{v}}^{(2)\to (1)}=\lambda \left(s\right){\mathbf{f}}^{(2)},$$

$${\mathbf{f}}^{(1)}+\lambda (s){\mathbf{f}}^{(2)}=-{\mathbf{J}}^{(1)}\cdot \left[\frac{\partial {\mathbf{r}}^{(1)}}{\partial t}-{\mathbf{U}}^{\infty }({\mathbf{r}}^{(1)})\right],$$

$${\mathbf{f}}^{(2)}+\lambda (s){\mathbf{f}}^{(1)}=-{\mathbf{J}}^{(2)}\cdot \left[\frac{\partial {\mathbf{r}}^{(2)}}{\partial t}-{\mathbf{U}}^{\infty }({\mathbf{r}}^{(2)})\right].$$

$${\mathbf{t}}^{(1)}=\cos \left(\theta -\varphi /2\right){\mathbf{e}}_x+\sin \left(\theta -\varphi /2\right){\mathbf{e}}_y,$$

$${\mathbf{t}}^{(2)}=\cos \left(\theta +\varphi /2\right){\mathbf{e}}_x+\sin \left(\theta +\varphi /2\right){\mathbf{e}}_y,$$

$${\mathbf{f}}^{(k)}(s)=-\frac{1}{1-{\lambda }^2(s)}\left({\mathbf{J}}^{(k)}\cdot \left[\frac{\partial {\mathbf{r}}^{(k)}}{\partial t}-{\mathbf{U}}^{\infty }({\mathbf{r}}^{(k)})\right]-\lambda (s){\mathbf{J}}^{(j)}\cdot \left[\frac{\partial {\mathbf{r}}^{(j)}}{\partial t}-{\mathbf{U}}^{\infty }({\mathbf{r}}^{(j)})\right]\right),k\ne j.$$

$${\mathbf{f}}^{(1)}+\lambda (s){\mathbf{f}}^{(2)}+\Lambda (s){\mathbf{f}}^{(3)}=-{\mathbf{J}}^{(1)}\cdot \left[\frac{\partial {\mathbf{r}}^{(1)}}{\partial t}-{\mathbf{U}}^{\infty }({\mathbf{r}}^{(1)})\right],$$

$${\mathbf{f}}^{(2)}+\lambda (s)({\mathbf{f}}^{(1)}+{\mathbf{f}}^{(3)})=-{\mathbf{J}}^{(2)}\cdot \left[\frac{\partial {\mathbf{r}}^{(2)}}{\partial t}-{\mathbf{U}}^{\infty }({\mathbf{r}}^{(2)})\right],$$

$${\mathbf{f}}^{(3)}+\lambda (s){\mathbf{f}}^{(2)}+\Lambda (s){\mathbf{f}}^{(1)}=-{\mathbf{J}}^{(3)}\cdot \left[\frac{\partial {\mathbf{r}}^{(3)}}{\partial t}-{\mathbf{U}}^{\infty }({\mathbf{r}}^{(3)})\right],$$

$${\mathbf{f}}^{(1)}$$

$${\mathbf{f}}^{(2)}$$

$${\mathbf{f}}^{(3)}$$

$$w_{k-1}=\exp \left[i\left(k-1\right)\varphi \right],1\le k\le N.$$

$$\sum _{k=1}^N\cos m{\theta }^{(k)}=\mathfrak{R}\left[\exp \left(im\theta \right)\sum _{k=1}^N{w}_{k-1}^m\right].$$

$$\sum _{k=1}^N{w}_{k-1}^m=\{\begin{array}{ll}0,& \text{if }m\text{ is not a multiple of }N,\\ N,& \text{if }m\text{ is a multiple of }N,\end{array}$$

Full text

1 Introduction

In microfluidic devices, manipulation of the flow and the suspended phases such as cells, droplets/bubbles, macromolecules (e.g. DNAs), etc. is commonly needed. The flow manipulation includes mixing, pumping, valving, sensing and related operations. In order to achieve these functions, different strategies have been developed. One of the most intuitive approaches is to introduce into the microfluidic device a rotary element, whose rotation is achieved by applying an external electric (Bart et al., 1992) or magnetic field (Ahn & Allen, 1995; Döpper et al., 1997; Ryu et al., 2004; Agarwal et al., 2005; van den Beld et al., 2015). On the other hand, such elements are also able to rotate passively without resorting to any external fields, but propelled by the underlying flow alone if they are placed asymmetrically with respect to the flow. This approach has been demonstrated by  Zaki et al. (1994) and  Day & Stone (2000) where the latter was inspired by the experimental work of rotating an asymmetrically placed cylinder to pump fluid in a duct (Sen et al., 1996). More recently,  Moon et al. (2015) has succeeded to drive rotary microgears by the underlying flow in a microfluidic channel. Each of their microgears consisted of eight paddles equally spaced in angle, and was fabricated and installed via in situ polymerisation based on flow lithography. The authors also showed that a pair of microgears was able to transmit the hydrodynamic torque from one gear to the other. Likewise, a similar flow-driven wheel was implemented by  Attia (2008) (PhD thesis in French) as a flow sensor to measure the flow speed in a microfluidic channel.

Motivated by such microfluidic experiments, we hereby carry out a theoretical and numerical study on the low-Reynolds-number hydrodynamics of such flow-driven rotary devices resembling micro-scale “watermills”. We aim to provide design principles for their applications in microsystems. After presenting the problem setup in Sec. 2, we describe in Sec. 3 the methodologies including the classical resistive force (RF) theory and numerical methods. The results obtained are compared in Sec. 4, which identifies the role of hydrodynamic interactions absent from the classical RF theory. Therefore, an improved RF theory taking hydrodynamic interactions into account is developed in Sec. 5.1 based on the recent theoretical work of  Man et al. (2016). We use the improved theory to solve the resistance and mobility problems of a simplified version of the watermill, and compare the theoretical predictions with the numerical results. Finally, we conclude and discuss our results in Sec. 6.

2 Problem setup

We consider a watermill-like rotary device consisting of $N\geq 1$ cylindrical rods equally distributed in a plane in the azimuthal direction (see figure 1). The angle between two successive rods is $\phi=2\pi/N$. All of the rods are jointed on a common end on the rotation axis of the device that is along the $z$ direction. Hence the device rotates in the horizontal $xy$ plane. The length of a rod is $L$, with its circular cross section of radius $a$, and the slenderness of the rod is defined as $\epsilon=a/L\ll 1$. We place the watermill in an unbounded Poiseuille flow $\mathbf{U}^{\infty}=U_{0}(1-x^{2}/R^{2})\mathbf{e}_{y}$, with the position $\mathbf{x}_{\mathrm{joint}}$ of the watermill’s joint (rotation axis) away from the centre ($x=0$) of the flow domain by distance $b\geq 0$, i.e., $\mathbf{x}_{\mathrm{joint}}=b\mathbf{e}_{x}$. Two nondimensional parameters

[TABLE]

are introduced to indicate the off-centre displacement of the watermill and the characteristic width of the channel flow, respectively. The orientation of the watermill is indicated by the angle $\theta$ between the $1$-st rod (arbitrarily labelled without losing generality) and the $x$-axis. The rotational velocity $\omega_{z}=\mathrm{d}\theta/\mathrm{d}t$ of the watermill depends on its orientation $\theta\in[0,\phi]$.

The dynamic viscosity of the fluid is $\mu$. We choose $U_{0}$, $L$, $L/U_{0}$, $\mu U_{0}/L$, $\mu U_{0}L$ and $\mu U_{0}L^{2}$ as, respectively the characteristic velocity, length, time, stress, force and torque. Nondimensional quantities are denoted by $\;\bar{}\;$ from hereafter. We fix $\epsilon=0.02$ in this study.

3 Methodologies

We carry out our study in the low-Reynolds-number flow regime and thus solve the Stokes equations. By employing the classical resistive force (RF) theory, we calculate the rotational velocities of a freely rotating watermill consisting of equally spaced rods. In Sec. 4, the results are compared with those computed by direct numerical simulations of the Stokes equations. The comparisons indicate that inter-rod hydrodynamic interactions neglected by the RF theory play an important role in certain configurations. This feature thus motivates us to conduct a theoretical study adopting the recently developed RF-based mathematical framework of  Man et al. (2016) that accounts for hydrodynamic interactions; this “RF-HI” theory will be described in Sec. 5. The classical RF theory and numerical methods are documented in Sec. 3.1 and  3.2, respectively.

3.1 Classical resistive force theory

We define an arclength $s\in[0,L]$ on each rod, with $s=0$ and $s=L$ corresponding to the joint and free end, respectively. The RF theory dictates that the hydrodynamic force per unit length $\mathbf{f}^{(k)}$ exerted by the fluid on the $k$-th rod is a function of the arclength $s$ according to

[TABLE]

where ${\mathbf{r}}(s)$ and $\mathbf{t}(s)=\partial{\mathbf{r}}/\partial s$ denote, respectively, the local coordinates and tangent of the rod’s centreline; $\mathbf{J}(\mathbf{t})=\xi_{\perp}\mathbf{I}+\left(\xi_{\parallel}-\xi_{\perp}\right)\mathbf{t}\mathbf{t}$, where $\xi_{\perp}$ and $\xi_{\parallel}$ are the drag coefficients for the motion of rod in the directions perpendicular and parallel to $\mathbf{t}$. Note that $\xi_{\perp}=2\xi_{\parallel}\approx 4\pi\mu/\ln{\epsilon^{-1}}$ (Lighthill, 1975).

The hydrodynamic torque exerted on the $k$-th rod about the joint is $\mathbf{T}^{(k)}=\int_{0}^{L}\left({\mathbf{r}}^{(k)}-\mathbf{x}_{\mathrm{joint}}\right)\times\mathbf{f}^{(k)}\mathrm{d}s$. Since the watermill only rotates in the $xy$ plane, we only consider the $z$-component of the torque on the $k$-th rod, whose nondimensional value is

[TABLE]

where $\bar{\omega}_{z}$ is the nondimensional rotational velocity of the watermill; $\theta^{(k)}$ denotes the angle between the $k$-th rod and the $x$-axis and $\theta^{(1)}=\theta$. The total torque on the rotary device, $\bar{T}_{z}=\sum_{k=1}^{N}\bar{T}_{z}^{(k)}$, from which we find

[TABLE]

Since the rods are equally spaced on a circle, we use the properties of roots of unity (detailed in appendix A) to obtain that

[TABLE]

[TABLE]

and

[TABLE]

so that the total torque $\bar{T}_{z}$ can be written as

[TABLE]

The rotational velocity of the freely rotating watermill can be obtained by applying the torque-free condition $\bar{T}_{z}=0$ and we find

[TABLE]

3.2 Numerical methods

To determine the rotational velocity of the watermill numerically, we carry out three-dimensional direct numerical simulations (DNS) based on a commercial finite-element method (FEM) solver COMSOL. We have experience in performing COMSOL simulations for viscoelastic flows, e.g. Pak et al. (2012) where the propulsion of two touching rotating spheres in viscoelastic fluids was investigated and the numerical results were in excellent agreement with the asymptotic analysis in the small Deborah number regime. It is worth noting that prior studies of interacting slender bodies in viscous flows have adopted other numerical implementations (Yamamoto & Matsuoka, 1995; Ross & Klingenberg, 1997; Saintillan & Shelley, 2007; Nazockdast et al., 2017).

Since we assume that the Reynolds number $Re=\rho_{f}U^{\infty}L/\mu$ is small ($\rho_{f}$ denotes the fluid density), inertia effects are negligible, and we solve the nondimensional steady Stokes equations

[TABLE]

For a channel flow, we impose a Dirichlet boundary condition (BC) with the Poiseuille flow profile in the inlet of the domain and a constant pressure BC at the outlet. Utilising the mirror symmetry of the setup, we only need to consider the upper half ($z\geq 0$) of the domain by applying a symmetry BC at the $z=0$ plane. We adopt the same symmetry BC at the $z=H_{z}/2$ plane, which effectively corresponds to an array of watermills equally spaced along the $z$ direction by distance $H_{z}$. No-slip BCs are specified on the two lateral walls at $x=\pm R$. Since both the RF and RF-HI theories are derived for unbounded flows, i.e. without accounting for wall effects, the boundedness of the computational domain needs to be considered carefully for a reasonable comparison of the theoretical and numerical results. We choose the length $H_{x}$ of the domain equal to $100L$. To mitigate the confinement effects of the lateral walls, the characteristic width of the channel flow $\rho=5$ is used in most of our cases. We have varied the distance $H_{z}$ in the range $\in[0.25,4]L$ and find $H_{z}=2L$ is sufficiently large to guarantee that the hydrodynamic interaction between the watermill and its mirror image about the symmetry BC at $z=H_{z}/2$ is negligible.

The BC imposed at position $\mathbf{x}_{w}$ on the surface $S_{w}$ of the watermill is

[TABLE]

where $\bar{\mathbf{x}}_{\mathrm{joint}}=\beta\mathbf{e}_{x}$ denotes the nondimensional position of the joint of the watermill. Note that the rotational velocity $\bar{\omega}_{z}$ is an unknown and it is solved together with the flow field $\left(\bar{p},\bar{\mathbf{u}}\right)$ by incorporating the constraint of zero hydrodynamic torque on the watermill

[TABLE]

where $\mathbf{n}_{w}$ denotes the unit normal vector on $S_{w}$.

The numerical setup is validated for a resistance and a mobility problem. For the resistance problem, we consider a cylindrical rod rotating at a constant velocity $\omega\mathbf{e}_{z}$ about one of its ends, with its revolution axis on the $z=0$ plane. The hydrodynamic torque $T_{z}/\left(\mu\omega L^{3}\right)$ calculated numerically agrees well with the RF predictions for varying slenderness $\epsilon$, as shown in table 3.2. For the mobility problem, we consider the free rotation of a spheroid in a shear flow, $\mathbf{U}^{\infty}=-\dot{\gamma}y\mathbf{e}_{x}$, where the revolution axis of the spheroid is on the shear plane ($z=0$). We define $2\kappa$ as the length of the axis of revolution scaled by the radius of a sphere with the same volume, and $\alpha$ as the angle between the revolution axis and the $x$-axis. We compute the instantaneous rotational velocity $\Omega=\Omega_{z}\mathbf{e}_{z}$ of the spheroid as a function of $\kappa$ and $\alpha$, and validate the results against the analytical theory (Jeffery, 1922), $\Omega_{z}/\dot{\gamma}=0.5\left[1-(\kappa^{3}-1)/(\kappa^{3}+1)\cos{2\alpha}\right]$. The comparison shown in table 3.2 shows a maximum discrepancy of below $1\%$.