Microfluidics involve the manipulation of samples and
reagents by utilizing extremely low volumes of fluids, down to the range of
femtolitres. A major application of microfluidics lab-on-chip devices involve
detecting, switching, and separating micro-entities - such as blood cells - in
micro-scaled channels^1,
2. These platforms have rekindled the interests of researchers to fully
understand the phenomenon of motion of suspended particles in a wall-bounded
environment. Today, a wide variety of Computational Fluid
Dynamics (CFD) packages are available to provide a solution for
non-linear continuum equations, reduce
the overall cost, and integrate complex components and processes in modeling
and optimizing lab-on-chip platforms³.
However, these CFD packages fail to capture a variety
of interesting phenomena, e.g., thermal fluctuations, particle-particle
interactions, and large-scale deformation of bioparticles in these complex
fluids at the mesoscale. Moreover, at such small scales, the discrete
nature of molecules becomes significant and cannot be ignored. Hence, we have
employed a particle-based technique called Dissipative Particle Dynamics
(DPD)

DPD is an off-lattice particle-based
method that relies on coarse-graining of molecules. Each DPD particle, modeled
by lumping a large number of molecules, moves together with other DPD particles
in a lagrangian fashion, interacting with them via soft repulsive
potentials. The inclusion of soft repulsive potential enables DPD to simulate
much larger length and time scales compared to MD. The DPD technique is
Galilean invariant, preserves both mass and momentum, and is capable of
producing correct hydrodynamics and colloidal behavior of the system⁴.
Presently, DPD has been successfully utilized to simulate the colloidal
suspensions, red blood cells in
capillaries⁵ and
microvessel stenosis⁶,
thrombosis formation⁷,
polymers, DNA, and multiphase flows^{8, 9}.

Newton’s second law governs the dynamical response
of the DPD particles. So, for a particle i:

Here,  r_i is the position
vector of the i^th DPD particle, ν_i represents
its velocity, is an external force acting on
the particle, and  is the internal force on
particle i due to its interaction with other DPD particles.
Specifically, three types of internal forces, effective within the cut-off
radius (r_c) exist in a DPD model: a soft, purely-repulsive conservative
force, a dissipative force to describe the viscosity of the DPD system, and a
random force to model
Brownian motion.  We have developed
an in-house FORTRAN code to model 2-D fluid flow followed by modeling
flow past a round obstacle in a microchannel. For Poiseuille flow, a
rectangular domain L×M is defined to represent the channel and the DPD fluid
particles with zero initial velocities are arranged randomly inside the domain.
Two layers of frozen wall particles are defined outside the fluid domain to
represent stationary microchannel walls. To model the flow around the
stationary obstacle, the 2-D round obstacle is modeled as a collection of
frozen DPD particles to form the circumference of obstacles. Each
circumferential particle is connected to its neighbors with a linear elastic
spring. The mass of both the fluid and obstacle particles are kept the same in
the simulation. To produce parabolic flow in the channel, a constant external
force is applied to all the fluid particles, and the modified velocity-verlet
integration algorithm introduced by Groot and Warren is used to solve their
equation of motions. To reduce the number of computations, a cell division and
link-list approach is utilized and only the interactions lying within the
cut-off radius are calculated.

A periodic boundary
condition is assigned to the left and right vertical boundaries to model an
infinitely long microchannel while the bounce back technique is utilized to
ensure no slip condition around the stationary channel walls and the obstacle
circumference. Although the fluid particles are bound to remain in the fluid
domain only; however, the soft potentials of the DPD technique allow the particles
to permeate the wall domain. If this happens, the leaving particle is
reinjected back to the fluid domain with reversed x- and y -velocity
components about the normal vector of channel wall or the obstacle boundary.
The DPD computational domain and the boundary conditions are shown in Fig. 1.

FIGURE 1. The
computational domain for flow past a stationary obstacle.

In the current work,
numerous parameters are changed to achieve a very small Re (= ρυD_H/µ) for the DPD
fluid flow. The cut-off radius was increased and the exponent of dissipative
weight function (s) is decreased to increase the dynamic viscosity of
the DPD fluid. The number density of the DPD fluid is kept 4. In order to
reduce these thermal fluctuations at low Re, the value of k_BT
value is reduced to recover a proper flow profile of the DPD fluid.

Figure 2 shows the velocity profile of the DPD fluid
at Re =
2.88, where Re is based on the channel height. The resulting
average velocity of the flow is 1.6975. The DPD result is compared with the analytical
result and a good agreement can be seen.

FIGURE 2. The (non-dimensional) velocity profile of
the DPD fluid at Re = 2.88.

The results of the fluid flow field around the 2-D
stationary round obstacle obtained from the current DPD model at Re= 3.2
is shown in Fig. 3. Here, Re is based on the diameter of the obstacle.
The blockage ratio – the ratio of the diameter of the obstacle to channel
height – is kept at 0.33 in
the current study.

FIGURE 3. The
streamwise velocity component (ν_x) past a frozen obstacle
(diameter = 5.28) in the DPD domain after 10⁴ time steps.

Finally, drag coefficient obtained from the drag force,
, on the stationary obstacle is obtained in DPD for Re
in the ranges of (~3 – 10). The drag force per unit length (F_d)
on the obstacle is calculated via vector summation of all the forces acting on
the circumferential beads of the obstacle, A is the cross-sectional area
perpendicular to the fluid flow, ρ is the fluid density, and  is the flow velocity relative to the object. To
verify the current DPD code, the DPD results are compared with already
published results¹⁰
and a close match is observed.

FIGURE 4.  Comparison of
the drag coefficient over a wide range of Reynolds number for the blockage
ratio of 0.33.

References

1.            W. Waheed, A. Alazzam, E. Abu-Nada,
S. Khashan and M. Abutayeh, Journal of Electrostatics 94, 1-7 (2018).

2.            A. Alazzam, B. Mathew and S.
Khashan, in Advanced Mechatronics and MEMS Devices II (Springer, 2017),
pp. 253-282.

3.            E. Fırat, presented at the AIP
Conference Proceedings, 2018 (unpublished).

4.            P. Hoogerbrugge and J. Koelman, EPL
(Europhysics Letters) 19 (3), 155 (1992).

5.            I. V. Pivkin and G. E. Karniadakis,
Physical review letters 101 (11), 118105 (2008).

6.            L. Xiao, S. Chen, C. Lin and Y. Liu,
Molecular & cellular biomechanics: MCB 11 (1), 67-85 (2014).

7.            N. Filipovic, M. Kojic and A. Tsuda,
Philosophical Transactions of the Royal Society of London A: Mathematical,
Physical and Engineering Sciences 366 (1879), 3265-3279 (2008).

8.            L. Yang and H. Yin, Physical Review
E 90 (3), 033311 (2014).

9.            P. De Palma, P. Valentini and M.
Napolitano, Physics of Fluids 18 (2), 027103 (2006).

10.          H. Reddy and J. Abraham, Physics of
Fluids 21 (5), 053303 (2009).