Mesoscale particle-based simulations of flow in expansion-contraction microchannels at low Reynolds number

TL;DR

This study compares analytical and mesoscale simulation methods for low Reynolds number flow in sinusoidal microchannels, finding good agreement and highlighting the strengths and limitations of each approach.

Contribution

It extends prior analytical solutions to tenth order and evaluates mesoscale particle-based methods against these solutions for complex microchannel flows.

Findings

DPD and MPCD simulations agree well with series solutions for flow rates and velocities.

Mesoscale models show slight wall slip, overpredicting flow parameters.

Series solutions fail for short channels and large amplitudes, but mesoscale methods remain effective.

Abstract

We computationally study the flow of Newtonian fluids through sinusoidal expansion-contraction microchannels at low Reynolds number. We first use a perturbation method to analytically derive series solutions for the stream function and volumetric flow rate that extend prior work [P.K. Kitanidis and B.B. Dykaar, Transport in Porous Media 26, 89-98 (1997)] up to tenth order. We then employ two particle-based mesoscale methods, dissipative particle dynamics (DPD) and multiparticle collision dynamics (MPCD), to simulate the same flows. We find that the fluid velocity at the expansion and contraction points as well as the volumetric flow rate are in good agreement between DPD, MPCD, and the fourth-order series solution for a wide range of microchannel geometries. The mesoscale fluid models exhibit some slip at the walls, leading to a small but consistent overprediction of the velocity and volumetric flow rate. The series solution fails for short microchannel lengths and large amplitudes; we identify lengths and amplitudes for which it converges to a given order. Overall, we find that DPD and MPCD are convenient and reasonably accurate methods, particularly for microchannel geometries where the series solution fails or is cumbersome to implement.

Equations

$$H(x)=\frac{W}{2}+A\cos \left(\frac{2\pi x}{L}\right),$$

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

$$\nabla \cdot \mathbf{\sigma }$$

$${\int }_S\mathbf{n}\cdot \mathbf{\sigma }\cdot \mathbf{u}d\mathbf{x}={\int }_V\mu \Phi d\mathbf{x},$$

$$Q={\int }_{A_{\mathrm{c}}}\mathbf{n}\cdot \mathbf{u}d\mathbf{x}$$

$$U_{\parallel }=\frac{W^2}{12\mu }\frac{\Delta p}{L}$$

$$\frac{{\partial }^4\hat{\psi }}{\partial {\hat{x}}^4}+2\frac{{\partial }^4\hat{\psi }}{\partial {\hat{x}}^2\partial {\hat{y}}^2}+\frac{{\partial }^4\hat{\psi }}{\partial {\hat{y}}^4}=0,$$

$$\hat{\Phi }=4{\left(\frac{{\partial }^2\hat{\psi }}{\partial \hat{x}\partial \hat{y}}\right)}^2+{\left(\frac{{\partial }^2\hat{\psi }}{\partial {\hat{y}}^2}-\frac{{\partial }^2\hat{\psi }}{\partial {\hat{x}}^2}\right)}^2.$$

$$\hat{H}(\eta )=\frac{1}{2}+\hat{A}\cos (2\pi \eta ),$$

$$\hat{Q}=\frac{1}{12}{\int }_{-1/2}^{1/2}{\int }_{-1}^1\hat{\Phi }\hat{H}d\xi d\eta .$$

$$\frac{\partial \hat{\psi }}{\partial \hat{x}}$$

$$\frac{\partial \hat{\psi }}{\partial \hat{y}}$$

$$\hat{\psi }(\eta ,\xi )=\sum _{n=0}^{\infty }{\hat{\psi }}_n(\eta ,\xi ){ϵ}^n,$$

$$\hat{Q}=\sum _{n=0}^{\infty }{\hat{Q}}_n{ϵ}^n$$

$${\hat{\psi }}_0$$

$${\hat{\psi }}_2$$

$${\hat{\psi }}_4$$

$${\hat{f}}_0$$

$${\hat{f}}_2$$

$${\hat{f}}_4$$

$$\displaystyle+\hat{H}^{2}\Big{(}\left(19-15\xi^{2}\right)\hat{H}\hat{H}_{4}+90\left(2\xi^{2}-3\right)\hat{H}_{2}^{2}\Big{)}$$

$$+8\left(30{\xi }^2-31\right){\hat{H}}^2{\hat{H}}_3{\hat{H}}_1$$

$$\displaystyle-36\left(50\xi^{2}-19\right)\hat{H}\hat{H}_{1}^{2}\hat{H}_{2}\Big{]},$$

$${\hat{Q}}_0$$

$${\hat{Q}}_2$$

$${\hat{Q}}_4$$

$${\hat{I}}_0$$

$${\hat{I}}_2$$

$${\hat{I}}_4$$

$$\displaystyle\quad+56\hat{H}^{2}\hat{H}_{3}\hat{H}_{1}-306\hat{H}\hat{H}_{1}^{2}\hat{H}_{2}\big{)}\hat{H}^{-3}{\operatorname{d}\!{\eta}}.$$

$$Q=\int Q(x)f^{(1)}(x)dx=\int A_{\mathrm{c}}(x)U(x)f^{(1)}(x)dx,$$

$$U(x)=\frac{\int v_x{f}^{(1)}(x,y,z,v_x)dydzd{v}_x}{\int f^{(1)}(x,y,z,v_x)dydzd{v}_x}.$$

$$f^{(1)}(x)=\int f^{(1)}(x,y,z,v_x)dydzd{v}_x$$

$$Q=⟨A_{\mathrm{c}}(x)v_x⟩.$$

Full text

I Introduction

Flow through expansion–contraction microchannels is used in a variety of technologies ranging from microfluidic mixers[1] to particle separators [2, 3, 4]. For example, expansion–contraction microchannels have been used to study DNA translocation [5, 6, 7] and the filtration of polymers with different architectures [8]. The curvature and constriction of expansion–contraction microchannels can trigger complex transport phenomena; for example, pulsatile forcing of a Newtonian fluid in a sinusoidal channel enhances mass transfer [9]. Expansion–contraction microchannels also serve as platforms for studying transport in porous media[10, 11, 12], especially the rich behaviors of non-Newtonian fluids in these geometries [13] that are relevant for processes such as enhanced oil recovery [14, 15, 16].

Both numerical and analytical methods have been applied to study flow in expansion–contraction geometries. Boundary integral methods were used to identify flow-reversal criteria and quantify geometric effects for Newtonian fluids in periodic constriction channels [17]. Pressure drop corrections and criteria for recirculation based on the amplitude-to-length ratio were derived with a similar method [18]. Lattice Boltzman simulations also have been applied to simulate Newtonian fluids in expansion–contraction geometries [19, 20]. Analytical perturbation methods were used to model Newtonian fluids in two-dimensional periodic pores [21] and three-dimensional wavy walls [22]. Non-Newtonian fluids are more challenging to describe with analytical approaches, but Boyko and Stone recently analyzed the behavior of Oldroyd-B type models in contractions using lubrication theory [23, 24, 25]. They also provide an excellent summary of prior numerical works on flow of non-Newtonian fluids in non-uniform geometries [23].

Mesoscale particle-based simulation methods such as dissipative particle dynamics (DPD) [26] and multiparticle collision dynamics (MPCD) [27] are intriguing alternatives for simulating microchannel flows. DPD represents a fluid as soft particles that interact through pairwise conservative, dissipative, and random forces, while MPCD represents a fluid as point-like particles that periodically participate in stochastic momentum-exchanging collisions. Both can efficiently model not only Newtonian fluids but also non-Newtonian fluids, e.g., by explicitly embedding coarse-grained polymers into a solvent. Such mesoscale representations may be critical for accurately describing transport phenomena for non-Newtonian fluids [28]. DPD has been used to study the trapping of fluid particles [29] and transport of macromolecules[30] in grooved nanochannels, flow of DNA in periodic contraction–diffusion channels[31], and polymer translocation dynamics [32, 33]; while MPCD has been used to investigate flow of colloidal particles in porous environments [34], translocation of knotted polymers [35], and filtration of ring polymers in expansion–contraction microchannels [36].

However, despite their popularity, to our knowledge there is no systematic study on whether DPD and MPCD accurately describe even the flow of Newtonian fluids in expansion–contraction microchannels. Establishing accuracy for this baseline case is important for using mesoscale methods to reliably simulate more complex phenomena, such as those for non-Newtonian fluids. In this work, we interrogate the ability of DPD and MPCD to model the flow of Newtonian fluids through sinusoidal expansion–contraction microchannels at low Reynolds number. The mesoscale particle-based simulations are compared to an approximate analytical theory. Overall, we find good agreement between both particle-based simulation methods and the theory for microchannel geometries for which the theory is expected to be accurate. Importantly, though, the particle-based methods can also be used to simulate flows in microchannels for which the theory fails.

The rest of this article is organized as follows. Section II presents the sinusoidal expansion–contraction microchannel geometry that we studied, along with a theoretical solution based on a perturbation method, while Section III describes the particle-based simulation models and methods. The simulation results are critically compared to theoretical predictions in Section IV; we examine the fluid velocity and volumetric flow rate in different microchannels, and we discuss shortcomings of both modeling approaches. Section V summarizes the main findings and provides an outlook.

II Theory

We considered an expansion–contraction microchannel that is effectively two-dimensional with a net flow in the $x$ direction and boundary walls in the $y$ direction (Figure 1). The half width between the walls $H$ varied sinusoidally as

[TABLE]

where $L$ is the microchannel length, $W$ is the average distance between the walls, and $A$ is the amplitude of the walls. The distance between the microchannel walls at the expansion and contraction points was hence $W\pm 2A$. The microchannel nominally also had depth $D$ in the $z$ direction, and there were periodic boundary conditions in the $x$ and $z$ directions. The domain for the fluid was taken to be $-L/2\leq x<L/2$, $-H(x)\leq y\leq H(x)$, and $-D/2\leq z<D/2$.

Kitanidis and Dykaar [21] studied the flow of a Newtonian fluid at low Reynolds numbers through a similar geometry and obtained a fourth-order series solution in the parameter $\epsilon=W/L$. Their approach is based on a stream function formulation of the Stokes equations, which they solved using a perturbation method [37]. The solution is expected to be accurate when the channel is long compared to the average wall separation, but they did not test its range of validity. Their derivation also has some notational inconsistencies and missing details. Here, we clarify how to solve for the flow in this microchannel geometry using their approach, and we extend their solution to tenth order for reasons that will be discussed in Section IV.

We assumed a steady incompressible flow of a Newtonian fluid with viscosity $\mu$ at low Reynolds number, so the fluid velocity $\mathbf{u}$ was governed by the Stokes equations [38],

[TABLE]

where $\boldsymbol{\sigma}=-p\boldsymbol{\delta}+\boldsymbol{\tau}$ is the stress tensor, $p$ is the pressure, $\boldsymbol{\delta}$ is the identity tensor, and $\boldsymbol{\tau}=\mu[\nabla\mathbf{u}+(\nabla\mathbf{u})^{\rm T}]$ is the deviatoric stress tensor. Written this way, the pressure may also include a conservative body force on the fluid. The boundary conditions were no slip and no penetration at the walls so $\mathbf{u}=\mathbf{0}$ at $y=\pm H$, periodicity in $x$ so $\mathbf{u}$ is equal at $x=\pm L/2$, and a constant pressure drop $\Delta p>0$ between $x=-L/2$ and $x=L/2$ to drive the flow in the positive $x$ direction. By symmetry, there should be no flow in the $z$ direction ($u_{z}=0$ everywhere), and $\mathbf{u}$ should not depend on $z$.

Kitadanis and Dyakaar [21] reformulated the pressure boundary condition as an auxiliary condition on the volumetric flow rate $Q$ by balancing the rate of work done on the fluid with the rate of viscous dissipation $\mu\Phi=\boldsymbol{\tau}:\nabla\mathbf{u}$ [39],

[TABLE]

where $S$ is the closed boundary of the microchannel volume $V$ with outer normal $\mathbf{n}$. We expected $\boldsymbol{\tau}$ to be periodic because it is a function of $\mathbf{u}$. Along with the boundary conditions and the integral definition of the volumetric flow

[TABLE]

through a cross-sectional area $A_{\rm c}$, the first integral was evaluated to be $\Delta pQ$. The second integral depends on the velocity $\mathbf{u}$ through $\Phi$.

The Stokes equations can be reformulated using a stream function $\psi(x,y)$ such that $\mathbf{u}=\nabla\times(\psi\mathbf{e}_{z})$, where $\mathbf{e}_{z}$ is the basis function in the $z$ direction. This stream function must satisfy the biharmonic equation $\nabla^{4}\psi=0$, and the boundary conditions on the velocity become boundary conditions on the corresponding first derivatives of $\psi$. For our microchannel geometry, $Q$ is also related to the stream function by Stokes’ theorem as $Q=[\psi(x,H)-\psi(x,-H)]D$. The stream function is only defined up to a constant so we arbitrarily chose $\psi(x,-H)=-Q/(2D)$, which implies $\psi(x,H)=Q/(2D)$.

To solve the biharmonic equation for $\psi$, we nondimensionalized the equations as follows: lengths were scaled by $W$, the stream function was scaled by $WU_{\parallel}$ where

[TABLE]

is the average velocity in a parallel-plate microchannel, and the volumetric flow rate was scaled by that of the parallel-plate microchannel, $Q_{\parallel}=WDU_{\parallel}$. All quantities in this section made dimensionless in this way will be denoted using a hat. Explicitly, the biharmonic equation became

[TABLE]

and we also had

[TABLE]

Next, we transformed the boundary conditions to use rescaled coordinates $\eta=\epsilon\hat{x}$ and $\xi=\hat{y}/\hat{H}$ having a rectangular domain $-1/2\leq\eta<1/2$ and $-1\leq\xi\leq 1$. In terms of these variables, the half width between the walls became

[TABLE]

and Eq. (4) became

[TABLE]

Some care must be taken in transforming the derivatives of $\hat{\psi}$ to use these coordinates. For example,

[TABLE]

where $\hat{H}_{i}$ designates the $i$-th derivative of $\hat{H}$ with respect $\eta$. Similar derivative transformations must be applied to all terms in Eqs. (7) and (8). The boundary conditions on $\hat{\psi}$ were $\hat{\psi}=\pm\hat{Q}/2$, $\partial\hat{\psi}/\partial\xi=0$, and $\partial\hat{\psi}/\partial\eta=0$ at $\xi=\pm 1$; and $\partial\hat{\psi}/\partial\xi$ and $\partial\hat{\psi}/\partial\eta$ being periodic at $\eta=\pm 1/2$.

An approximate solution to the biharmonic equation with these boundary conditions can be derived using a perturbation method with the small parameter $\epsilon$, i.e., microchannels that are long compared to their average width [21]. We sought solutions of the form

[TABLE]

where $\{\hat{\psi}_{n}\}$ are unknown functions to be determined. Given the form of Eq. (10), we also expected

[TABLE]

to have a corresponding series, where $\{\hat{Q}_{n}\}$ are unknown coefficients to be determined. The series for $\hat{\psi}$ was substituted in the biharmonic equation, then like-powers of $\epsilon$ were collected to obtain a set of partial differential equations for $\{\hat{\psi}_{n}\}$. These partial differential equations could be solved using only the boundary conditions $\partial\hat{\psi}_{n}/\partial\xi=0$ at $\xi=\pm 1$ and $\hat{\psi}_{n}=\pm\hat{Q}_{n}/2$ at $\xi=\pm 1$. The solutions for $\{\hat{\psi}_{n}\}$ contained the unknown coefficients $\{\hat{Q}_{n}\}$, which we subsequently determined by substituting $\{\hat{\psi}_{n}\}$ in Eq. (10) and solving a system of linear equations. We performed these calculations using Mathematica version 14.2.1.0, which enabled us to obtain a solution up through tenth order.

The solution contained only even powers of $\epsilon$. The terms through fourth order were

[TABLE]

with functions $\{\hat{f}_{n}(\eta,\xi)\}$,

[TABLE]

and coefficients

[TABLE]

along with the integrals

[TABLE]

We verified these results were equivalent to those of Kitanidis and Dyakaar with different notation[21]. We present complete details of the calculations through tenth order in the Supplementary Material.

III Simulation

We simulated the same microchannel geometry described in Section II using two mesoscale particle-based methods: DPD [40, 41] and MPCD [42]. To facilitate comparison of these methods, all simulation parameters and results will be presented using a consistent system of units having unit of mass $m$, unit of length $\ell$, and unit of energy $\varepsilon$. The unit of time is $\tau=\sqrt{m\ell^{2}/\varepsilon}$. All particles had mass $M=1.0\,m$, and all simulations were conducted at constant temperature $T=1.0\,\varepsilon/k_{\rm B}$, where $k_{\rm B}$ is the Boltzmann constant. Note that $\varepsilon$ used here and later in this article as the unit of energy should not be confused with $\epsilon=W/L$ in Section II.

Simulations were performed for microchannel lengths $L$ between $30\,\ell$ and $150\ell$. The mean width between walls $W$ and the depth $D$ were fixed at $30\,\ell$, and the wall amplitude $A$ was varied between $0\,\ell$ and $12\,\ell$. Periodic boundary conditions were applied in the $x$ and $z$ directions, and a bounce-back rule was used to model no-slip and no-penetration boundary conditions at the walls in the $y$ direction [43, 44]. Intersection of particles with the walls was computed numerically using a combination of Newton’s method and bisection search, and the velocities of particles that intersected with the wall were reflected. The orthorhombic simulation box containing the microchannel was padded by $4\,\ell$ in both the positive and negative $y$ direction at the expansion point to accommodate wall particles for DPD [45, 46] and virtual particles [43, 44] for MPCD (see below). All simulations were performed using HOOMD-blue [47, 48, 49] (version 5.3.0) extended with azplugins [50] (version 1.1.0). We present specific details of the two simulation methods next.

III.1 Dissipative particle dynamics

DPD uses particles that have three types of pairwise interactions: a conservative force, a dissipative force, and a random force [40, 41, 26]. In our simulations, all particles interacted through the standard DPD conservative repulsive force with interaction parameter $25\,\varepsilon/\ell$, but we used a modified weight function with exponent $1/2$ and friction coefficient $4.5\,m/\tau$ for the dissipative and random forces to produce a more liquid-like fluid [31, 51]. The pairwise distance cutoff for all forces was $1.0\,\ell$, and the equations of motion for the particles were integrated using a velocity Verlet scheme[52] with timestep $0.01\,\tau$.

For each microchannel geometry, particles were first randomly dispersed everywhere in the simulation box with no walls and three-dimensional periodic boundary conditions at a number density of $\rho_{\rm n}=3/\ell^{3}$. This bulk fluid was equilibrated for $10^{3}\tau$ then converted to represent a fluid confined between walls. Particles inside the microchannel continued to represent the fluid [blue particles in Fig. 1(a)], and their center-of-mass velocity was set to zero. Particles outside the microchannel but within a distance of $3\,\ell$ in the $y$ direction were converted to be wall particles [gray particles in Fig. 1(a)], their velocities were set to zero, and their motion was no longer integrated. All other particles were deleted. The presence of explicit wall particles helped reduce slip and undesirable density oscillations compared to using only bounce-back reflection [45, 53, 46, 51]. The fluid particles were then equilibrated for an additional $10^{3}\,\tau$ before starting a flow simulation.

III.2 Multiparticle collision dynamics

MPCD uses particles that do not interact through pairwise forces; instead, their positions and velocities evolve through alternating streaming and collision steps [42, 55, 27]. In our simulations, the particles moved ballistically for $0.1\,\tau$ during the streaming step. Their positions and velocities were updated using a velocity Verlet scheme [52], and bounce-back reflection was enforced with the walls [43, 44]. During the collision step, the particles were binned into cubic cells of edge length $1.0\,\ell$ then exchanged momentum with other particles in the same cell using the stochastic rotation dynamics scheme [42] with a fixed rotation angle of $130^{\circ}$, a rotation axis randomly drawn from the unit sphere for each cell, and a cell-level Maxwell–Boltzmann thermostat [56]. Following best practices, the cells were shifted in each Cartesian direction by a uniform random amount between $-\ell/2$ and $\ell/2$ before each collision [57]. Virtual particles were also added everywhere outside the walls using random sampling with rejection and participated in the collision [43, 44].

For each microchannel geometry, MPCD particles were randomly inserted between the walls at number density $\rho_{\rm n}=10/\ell^{3}$, and a flow simulation was immediately commenced.

III.3 Simulation protocol

To simulate flow with both methods, a constant force $F$ was applied to all fluid particles in the $x$ direction [58]. In selecting the value of $F$, we desired to have a small Reynolds number ${\rm Re}=\rho U_{\parallel}W/\mu$, where $\rho=M\rho_{\rm n}$ is the mass density of the fluid and $U_{\parallel}$ is calculated using Eq. (6) with $\Delta p/L=\rho_{\rm n}F$, for consistency with physical microchannel flows and the assumption underlying our theoretical analysis (Section II) yet a large enough velocity to be reliably measured. We used $F=10^{-3}\,\varepsilon/\ell$ in the DPD simulations and $F=2.5\times 10^{-3}\,\varepsilon/\ell$ in the MPCD simulations, which gave ${\rm Re}\approx 7.6$ and 7.5, respectively, based on the viscosities we calculated (Section IV). In each simulation, the flow was allowed to develop for $10^{4}\tau$, then the fluid velocity $\mathbf{u}$ and volumetric flow rate $Q$ were measured using the procedures described below during a production period of $10^{5}\tau$.

The fluid velocity at the expansion and contraction points $u_{x}(y)$ was measured using histograms sampled every $10\,\tau$. Fluid particles with an $x$ coordinate within a thickness $0.5\,\ell$ centered at the expansion were binned with respect to their $y$ coordinate between the walls using a bin size of $0.5\,\ell$, and the velocities of particles in each bin were averaged. The same procedure was used at the contraction, but due to the periodic boundary conditions, particles were binned in two separate half-thicknesses at each end of the simulation box and the resulting histograms were averaged together with equal weight. The two-dimensional velocity field $\mathbf{u}(x,y)$ was also calculated by binning the entire simulation box with respect to the $x$ and $y$ coordinates of the fluid particles using the same bin size and sampling frequency.

The volumetric flow rate $Q(x)$ at a given position $x$ can be calculated using Eq. (5) by integrating $\mathbf{u}(x,y)$ over the cross-sectional area $A_{\rm c}(x)=2H(x)D$. However, the MPCD fluid can exhibit significant density variations in response to pressure variations (Figure S2) because it has an ideal-gas equation of state [42, 59, 60, 61]. Variations in density should lead to variations in $Q$ because of mass conservation, so a volumetric flow rate averaged over the microchannel may be more appropriate than averaging only at a particular point (e.g., the expansion or contraction). Accordingly, we instead computed the volumetric flow rate using an ensemble average over the particles.

Based on Eq. (5), we defined the average volumetric flow rate as

[TABLE]

where $f^{(1)}(x)$ is the one-particle distribution function for finding a particle at position $x$[62], and $U(x)$ is the average $x$-component of the fluid velocity through the cross-sectional area $A_{\rm c}(x)$ at $x$. We note that there is not a unique definition of this average and others are possible; we used this one because it more heavily weights $x$ positions where fluid is more likely to be found. We then used the one-particle distribution function $f^{(1)}(x,y,z,v_{x})$ for finding a particle at position $(x,y,z)$ with $x$-velocity $v_{x}$ to express the average velocity as

[TABLE]

Using the fact that

[TABLE]

and simplifying gave the ensemble average

[TABLE]

We performed this average in our simulation by sampling the particle coordinates every $10^{2}\,\tau$.

IV Results and Discussion

We compared the fluid velocity and volumetric flow rates simulated using DPD and MPCD (Section III) with those calculated theoretically (Section II). To make this comparison, we required the viscosity $\mu$ of the DPD and MPCD fluids. Expressions for $\mu$ in terms of model parameters exist for both methods [63, 31]; the expression for MPCD is known to be quite accurate [64], but the expression for DPD is less accurate [51]. Accordingly, we determined $\mu$ for both methods using our simulations.

To use as much of our data as possible, we averaged the two-dimensional velocity in two parallel-plate microchannels ($A=0\,\ell$, $L=75\,\ell$ and $150\,\ell$) with respect to $x$ using equal weight to obtain $u_{x}(y)$ (Figure S3). We expected a parabolic flow $u_{x}(y)=(3U_{\parallel}/2)[1-(2y/W)^{2}]$ in this geometry for no-slip boundary conditions. However, to avoid ambiguities in determining the viscosity if there is wall slip, we used the velocity gradient ${\operatorname{d}\!{u_{x}}}/{\operatorname{d}\!{y}}=-12U_{\parallel}y/W^{2}$. We numerically differentiated the measured $u_{x}$ with respect to $y$; it was mostly linear, as expected, but there were small deviations near the walls (Figure S4). We fit ${\operatorname{d}\!{u_{x}}}/{\operatorname{d}\!{y}}$ for $|y|\leq 10\,\ell$ to a linear form, and we extracted viscosities of $1.63\,\varepsilon\tau/\ell^{3}$ for DPD and $8.68\,\varepsilon\tau/\ell^{3}$ for MPCD. The viscosity for DPD is in good agreement with previous measurements using reverse-nonequilibrium simulations [51] but is about 14% larger than the theoretical prediction of $1.43\,\varepsilon\tau/\ell^{3}$ [31], while the viscosity for MPCD is in excellent agreement with the theoretical expectation of $8.70\,\varepsilon\tau/\ell^{3}$ [63]. We used the values we measured to normalize the simulation results in all subsequent comparisons.

IV.1 Velocity

We first compared the simulated fluid velocities at the expansion and contraction points to the fourth-order theoretical predictions for microchannels with fixed length $L/W=5$ and amplitudes varying from $A/W=0$ to 0.4 (Figure 3). The theoretical predictions for the velocity were obtained by analytically differentiating the stream function. The maximum velocity at the expansion systematically decreased as $A$ increased; at the contraction, it initially increased before decreasing. In all cases, the DPD and MPCD simulations were in nearly quantitative agreement, demonstrating consistency of the methods, but both gave velocities that were slightly larger than the theoretical predictions (Figure 3 and S5). We expected the theoretical predictions to be accurate for these microchannel geometries because they were long compared to their width ($W/L=0.2$ is reasonably small). We also saw a comparable amount of disagreement between the simulations and theory even for the parallel-plate microchannel, for which the theoretical prediction is exact. We hence attribute the difference between the simulations and theory to the presence of a small amount of effective wall slip.

In support of this notion, we refit the velocities in the parallel-plate microchannels that we used to determine the viscosity, but we now allowed for a finite slip velocity at the walls (Figure S3). The fit slip velocity was 2.3% of $U_{\parallel}$ for DPD and 3.2% for MPCD. This amount of slip is consistent with deviations in the velocity for all microchannels in Figure 3. Some possible origins of the effective slip might be imperfect enforcement of no-slip boundary conditions using the bounce-back reflection scheme, as well as small deviations in fluid properties near the walls. However, the agreement between simulations and theory is overall very good, validating the implementation of both simulation models for the expansion–contraction microchannel as well as the theoretical calculations.

IV.2 Volumetric flow rate

Given the good agreement in the velocity between the simulations and theory for long microchannels ($L/W=5$), we next tested whether similarly good agreement could be obtained in shorter microchannels. We halved the length to $L/W=2.5$, and we conducted simulations for the same amplitudes as in Figure 3. In order to more easily visualize the results for multiple microchannels, we will now focus on the volumetric flow rate $Q$ rather than the velocity [Figure 4(a)]. For both lengths, $Q$ decreased as the amplitude increased due to the enhanced constriction of the flow, reaching $Q/Q_{\parallel}\approx 0.1$ for the most constricted channels. Both the DPD and MPCD simulations were again in good agreement with the fourth-order theoretical predictions, with the simulated values being between 2.5% (for $A/W=0$) and 8.5% (for $A/W=0.4$) larger than the theoretical predictions when $L/W=5$. This overprediction is consistent with the effective slip we noted in Figure 3.

To assess whether similarly good agreement could be achieved for other microchannel lengths, we fixed the amplitude at $A/W=0.2$ and 0.4, and we varied the length across the range $1\leq L/W\leq 5$ [Figure 4(b)]. For these geometries, $Q$ tended to decrease as $L$ decreased, qualitatively because the fluid encounters constrictions more frequently when $L$ is smaller. As before, the DPD and MPCD simulations produced almost identical results, and we found again that there was good agreement between the simulations and theoretical predictions when $L/W\gtrsim 2$. For shorter channels, however, the theoretical predictions began to deviate noticeably from the simulations. We suspected that this deviation was due to inaccuracy in the theoretical predictions, rather than the simulations, because the small parameter used in the perturbation method increases as $L/W$ decreases, potentially making the series expansion less accurate.

To test this idea, we calculated the theoretical solution for the volumetric flow rate using approximations from zeroth to tenth order as a function of length $L$ at fixed amplitude $A/W=0.2$, and we compared these predictions to the simulation results [Figure 5(a)]. Increasing the order of the approximation improved the match to the simulations at moderate lengths, but an apparent divergence emerged for small lengths; the predicted $Q$ seemed to alternate between increasingly positive and negative values for $L/W=1$. To illustrate this behavior more clearly, we also show the theoretical predictions for $Q$ as a function of the approximation order in short microchannels ($1\leq L/W\leq 2$) when $A/W=0.2$ [Figure 5(b)]. We found that $Q$ tended to converge when $L/W\geq 1.4$, was nearly oscillatory when $L/W=1.2$, and was clearly divergent when $L/W=1$.

To provide guidance on use of the theoretical solution, we identified which geometries could be modeled by what order series (Figure 6). As a simple test, we considered a series converged at order $n$ if its sum was positive for all orders of approximation and if all series terms of order $n$ and greater had absolute value less than 1% of the tenth-order sum; otherwise, the series was considered unconverged. As expected and consistent with Figure 5, the order necessary for convergence monotonically increased as the microchannel length decreased at fixed amplitude. Further, for a fixed length, microchannels with larger amplitudes typically required a higher-order series to achieve convergence. There is also a region of short lengths and large amplitudes (upper left corner) where the series did not converge even to tenth order.

V Conclusions

We investigated the flow of Newtonian fluids through sinusoidal expansion–contraction microchannels at low Reynolds numbers using complementary analytical and numerical methods. Extending prior analysis [21], we first derived series solutions for the stream function and volumetric flow rate in these geometries up to tenth order using a perturbation method. We then used two popular mesoscale particle-based fluid models, DPD and MPCD, to simulate the same flows numerically. Overall, we found that the particle-based simulations and the series solution were in good agreement with each other. The fluid velocity at the expansion and contraction points not only matched between DPD and MPCD, but also agreed with the fourth-order series solution. The volumetric flow rate also showed good agreement for a range of microchannel geometries. We did note a small overestimation of the flow rates in both particle-based methods, which we attributed to the presence of some effective wall slip in these methods.

For shorter microchannels, we found noticeable deviations between the theoretical solution and the particle-based simulations. We showed that these deviations were due to issues with convergence of the series solution. We also identified channel geometries for which the series solution for the volumetric flow rate converged for a given order of approximation. This analysis can be used as a guide for which level of approximation is needed for a specific channel geometry. Mesoscale particle-based methods like DPD or MPCD may be particularly useful for geometries where the series solution fails to converge or requires high-order terms that are cumbersome to implement. Additionally, both DPD and MPCD have straightforward ways [27, 65] to incorporate various solutes, such as polymers, to simulate complex fluids using the same approach. Such extensions may enable new computational studies of viscoelastic phenomena in expansion–contraction geometries such as elastic turbulence [13].

Supplementary Material

See the supplementary material for alternative normalization of velocities shown in Figure 3, density distribution and two-dimensional velocity for MPCD fluid, parallel-plate velocity and gradients used to extract viscosities, and complete details for series solutions for stream function and volumetric flow rate up to tenth order.

Conflicts of interest

The authors have no conflicts to disclose.

Data Availability

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

Acknowledgements.

Acknowledgment is made to the donors of the American Chemical Society Petroleum Research Fund for support of this research (Grant Nos. 65334-DNI7 to AS and 66616-DNI9 to MPH). This work used Delta at the National Center for Supercomputing Applications through allocation CHM240073 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program [66], which is supported by U.S. National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296.