1 Introduction

Arrays of micropillar are gaining popularity quickly due to their potential applications in a wide range of areas. Flow across arrays of pillars embedded inside microfluidic chip systems is one of the most active research fields and is of great importance for high-performance liquid chromatography^[1], thermal management^[2], dielectrophoresis^[3], and isolating diseased cells^[4]. These micropillars can increase mass and heat transfer coefficients, surface-to-volume ratio, surface chemical reactions, and heat conductivity with the compromise of a decrease in the effective cross section area^[5]. The decrease of the flow area causes an overall increase of the pressure drop. This competitive trend indicates that an optimum value can be found for various application designs. One of the important parameters for arrays of micropillar inside microfluidic chips is the volume fraction of micropillar or porosity (fluid volume fraction). Small porosity with large pillar radii or densely alignment can lead to a high driving capillary pressure and decrease the permeability. For micropillar arrays to be considered as an effective porous medium, it is very important to accurately predict the maximum mass flow rate pass arrays based on the micropillar array design.

Flow around regularly arranged micropillar has been the subject of numerous analytical and numerical works. Most of numerical works^[6-9] focus on the analytical solution of the flow resistance to the two-dimensional flow through a periodic array of circular cylinders. In experimental or analytical works, it is difficult to accurately predict the flow through the cylinder array when the cylinder array embedded inside microfluidic chip and channel boundary effects cannot be dismissed^[5]. As an effective alternate, numerical simulation of flows in microfluidic with pillar arrays attracted more and more researchers. Srivastava et al.^[10] and Ranjan et al.^[11] used the finite element method to simulate flow through a square pillar array, and they derived a function to describe the relationship between the permeability and the dimensionless geometric parameters. Ellero and Adams^[12] presented a numerical approach based on the smoothed particle hydrodynamics (SPH) to study a very viscous Newtonian liquid through a linear cylinder array embedded in a channel. Liu et al.^[13] presented a modified dissipative particle dynamics (DPD) model to study the unsaturated multiphase flow through porous media.

This paper presents a particle-based simulation approach, the DPD^[14], to model flow through periodic arrays of circular micropillar. The DPD is a recently developed mesoscale approach which can be used to simulate the hydrodynamic behavior of simple and complex fluids. As Hoogerbrugge and Koelman^[14] stated, the DPD method is a novel particle-based approach combining the good features of both the molecular dynamics (MD) and lattice gas automata (LGA) approaches, which is much more flexible than the LGA and much faster than the MD. In the DPD framework, each particle represents a cluster of atoms or molecules that interact via conservative (non-dissipative), dissipative, and random forces. The repulsive potential between the DPD particles is much softer than that between individual atoms or molecules, making larger time steps possible relative to the MD approaches. The longer time steps and the larger particle size can be used in the DPD, which makes the DPD much more practical than the MD to simulate complex hydrodynamics. Because the number of particles in the DPD system remains unchanged and each particle has the same mass, the DPD conserves the total mass of the system. The total momentum in the typical DPD simulation is conserved because the interactions between the particles are symmetrical. The Galilean invariant is observed, since the interactions between particles only depend on relative positions and velocities. The DPD is more flexible than the lattice Boltzmann method (LBM) since it does not suffer from the numerical instability. It is not as computationally efficient as the LBM simulations. It is more suitable to use the DPD for the simulation of complex fluids on interestingly physical time and length scales, including polymer suspensions^[18], colloids and gels^[20], and liquids with interfaces^[19].

In this work, we present a DPD model to simulate flow through periodic arrays of circular micropillar. The paper is organized as follows. In Sections 2 and 3, the DPD methodology and theory for circular micropillar arrays are described, respectively. In Section 4, the results of numerical simulations are provided and analyzed. The paper concludes in Section 5 with some remarks.

2 DPD method

In the DPD system, it is generally assumed that all particles have equal mass, and use this mass as a unit. The interactions between particles are assumed to be pairwise additive. The time evolution of DPD particles is governed by Newton’s equations of motion. Therefore, the governing equation for a simple DPD particle i can be expressed as

(1)

where r_i and v_i are the position and the velocity of a particle i, respectively, m_i is the mass of a DPD particle i, usually taken as unity, and f_i denotes the total force acting on a particle i. f_i^ext is the external force, such as the gravity. The inter-particle force acting on a particle i, f_i^int consists of three parts, namely, the conservative force F_ij^C, the dissipative force F_ij^D, and the random force F_ij^R,

(2)

where F_ij denotes the inter-particle force on the particle i by the particle j, which is equal to F_ji in the magnitude but opposite in the direction. The pairwise inter-particle interactions have a finite cut-off radius, r_c, which is generally taken as the length unit in the DPD systems.

The conservative force describes the thermodynamic properties of the fluid and can be derived from the interaction potential between the particles i and j as

(3)

where a_ij represents the repulsion amplitude between the particles i and j; r_ij=r_i -r_j represents the relative position between the particles i and j; r_ij=|r_ij|, ȓ_ij=r_ij/r_ij is the unit vector directed from the mass center of the particles j to i; w^C(r_ij) is the conservative weight function, and in the conventional DPD implementations, it takes a simple form as

(4)

Here, r_c is the cut-off distance. Equations (3) and (4) show that the conservative force is only repulsive over a finite distance r_c and acts in the direction of the ȓ_ij vector. The soft conservative force makes it possible to choose a large time step in a numerical method for integrating the equations of motion.

The dissipative force represents the viscous effects of the DPD system, while the random force is responsible for the eliminated degrees of freedom. Both dissipative and random forces are dependent on the relative positions of the particles. Moreover, the random force is related to relative velocities of the particles. The dissipative force and the random force are written as

(5)

(6)

where γ and σ are two coefficients which describe the strength of the dissipative and random forces, respectively. w^D (r_ij) and w^R(r_ij) are two weight functions which represent the variation of the friction coefficient and the noise amplitude with distance, respectively. v_ij (=v_i — v_j) is the relative velocity between the particles i and j. ξ_ij is a random variable with Gaussian statistics.

It is seen that the dissipative force acts to reduce the relative velocity of the particles, thus removing the kinetic energy and decreasing the temperature of the system, while the random force acts to heat up the system. Therefore, the dissipative and random forces act together to keep the system at a constant temperature with small fluctuations about the designed tem¬perature T. In order to correct isothermal balance, the coefficients (γ and σ) and the weight functions (w^D and w^R) of the random and dissipative forces must satisfy two constraints^[16],

(7)

(8)

where k_B is the Boltzmann constant, and T is the absolute temperature of the system. k_BT is usually taken as the energy unit, and then it can be used to express all the interaction energies of the system. There are different forms of w^D(r) and w^R(r). In a typical DPD formulation, w^D(r) and w^R(r) are

(9)

where r_d is the cut-off radius for the dissipative and random forces. The choice of r_d may be different from that of the conservative force r_c, and it influences the computation cost by affecting the number of pairs of interacting particles. s is the exponent of the weighting function, while different s can result in different dynamic behaviors of the DPD system. It is shown in Refs. [15] and [17] that s can influence the viscosity and the Schmidt number of the fluid by affecting the strength of dissipative force between particles. In a typical DPD simulation, the choices can be taken as r_c=r_d=1 and s=2.

3 Theoretical analysis

The correct simulation of flow around a single circular pillar is the base of modeling flow through multi pillars. To this end, we first use the DPD to simulate flow around a circular pillar. Then, the permeability is calculated for circular pillar arrays with different porosities. In this section, the theories for flow around single circular pillar and circular pillar arrays are presented, respectively, and further a comparison with simulation results is presented in next section.

3.1 Flow around single circular pillar

Before modeling the flow around circular pillar arrays, we first simulate the flow around a single stationary circular pillar between two infinite parallel plates. As shown in Fig. 1, the distance between two plates is h. The density and the viscosity of fluid are ρ and μ, respectively. The pillar is at the center between two plates, and its diameter is d.

The hydrodynamic drag force exerted on a rigid object suspended in a Newtonian fluid flow has been intensely researched over the past century. The drag force of an object depends on the properties and velocity of the fluid and the size and shape of the object. The drag force of a stationary object can be described as follows:

(10)

where v is the velocity of the fluid. A is the reference area of the object, and C_D is the drag coefficient. Consider that the diameter of the pillar is very small compared with the distance between two plates, namely, d < < h, making it possible to ignore the h effect. The experiment^[21] showed that the drag coefficient is a function of the Reynolds number for two-dimensional circular cylinders, and in small Reynolds numbers, the drag coefficient can be approximated as follows:

(11)

3.2 Flow around circular pillar arrays

In this paper, we are finally interested in calculating the permeability of micropillar arrays. Two types of periodic circular micropillar arrays are shown in Figs. 2 and 3. The pillar diameter is d, and the solid volume fraction of the arrays is ϕ=πd²/(4L²). Spielman and Goren^[22] obtained a function of permeability for fibrous arrays by simulation, which has been proven^[23-25]. The function gives relationships between the dimensionless permeability K/d2 (K is the permeability) and the solid volume fraction of the medium ϕ. For matrices with fibers aligned perpendicular to the flow, the form of the function is

(12)

where K_n(x) is the modified Bessel function of the second kind and order n with argument x, and ϕ < 0.75. From Eq. (12), it should be noted that the dimensionless permeability is an immutable value when the solid volume fraction is constant. In the present DPD simulation, the flow is driven by a body force g exerting on every DPD particle. Therefore, the permeability can be expressed as

(13)

4 Simulation and results

In Subsection 4.1, we use the DPD to simulate the flow around a single micropillar. The relationship between the drag coefficient and the Reynolds number Re of the experiment is compared with that of the DPD simulation. In Subsection 4.2, we use the DPD to simulate the flow around micropillar arrays, and the simulation permeability of circular micropillar arrays is compared with Eq. (12).

4.1 Flow around single micropillar

In this subsection, the DPD is used to simulate the flow around a single micropillar with parameters ρ=4, r=1, k_BT=1, γ=4.5, σ=3, a=18.75, and s=0.25. With these parametric values, the viscosity of DPD fluid is μ=2.40. As shown in Fig. 1, the height and length of the channel are h=30 and l=60, respectively. The diameter of the pillar is d=6. The non-slip boundary condition is implemented by frozen particles in the solid obstacle areas together with the Maxwellian reflection model when a fluid particle enters a thin reflecting boundary layer. This solid boundary treatment has been proven to be effective in preventing the unphysical particle penetration, and density oscillations can be controlled in a reasonably low level^[15-26]. In our simulation, different Reynolds numbers are simulated by adding different body forces g exerting on every DPD fluid particle. In the xz-plane, the computational domain is divided into 120 × 66 bins, and local data are collected in each bin. We can obtain all local flow properties by averaging the sampled data in each bin over 10⁵ time steps. Horizontal and vertical velocity fields are shown in Figs. 4 and 5. The results indicate that the horizontal velocity tends to zero near the channel walls and pillar. The vertical velocity is almost zero in the region far away from the cylinder, while it increases in the region above or below the cylinder. The number density and streamline of the flow field are shown in Figs. 6 and 7. The number density is almost uniform across the channel. The streamline of the flow field indicates that the velocity along the surface of the cylinder is parallel to the surface of the cylinder, which is consistent with the horizontal velocity field and the vertical velocity field.

The flow past a two-dimensional cylinder is one of the most studied problems in hydro¬dynamics. It is relevant to many engineering applications. Studies have shown that the flow pattern and the drag on a pillar under steady flow only related to the Reynolds number. In our simulation, the drag force is determined by computing the vector sum of the forces acting on all particles comprising the cylinder or micropillar, and its value is averaged over 10⁵ time steps. The comparison of the drag coefficient of circular cylinder is shown in Fig. 8. Figure 8 shows that the drag coefficient of circular cylinder agrees well with the experimental results for Re from 0 up to 100. The error between the numerical and experimental data increases with the increase of Re. When Re is increased to 100, the drag coefficient from the simulation is smaller than that from the experiment.

4.2 Flow around circular micropillar arrays

After successful simulation of the flow around a single micropillar, we use the DPD to simu¬late the flow around two types of circular micropillar arrays in this section. The computational domain is indicated in the dashed frame, as shown in Figs. 2 and 3. All DPD parameters are the same as those mentioned in the previous simulation. By changing the distance between two consecutive pillars with fixed diameter, we can obtain the differential permeability. In our simulation, different Reynolds numbers are simulated by adding different body forces g exerting on every DPD fluid particle. The flow fields including the horizontal and vertical velocity fields, the number density field, and the streamline of the flow of square and hexagonal arrays are shown in Figs. 9 and 10, respectively. The relationships between the horizonal body force and the mean horizonal velocities of square and hexagonal arrays are shown in Figs. 11 and 12, respectively. From Figs. 11 and 12, we can see that the mean horizonal velocities of both square and hexagonal arrays are proportional to the horizonal body force, which is consistent with the theoretical equation (13). In legends of Figs. 11 and 12, ϕ and ϕ' represent the actual value and the predicted value of solid volume fraction, respectively. ϕ is calculated from the initial computational setup, and ϕ' is estimated by putting permeability from the DPD simulation into Eq. (12). The predicted solid volume is quite close to the actual counterpart, especially the moderate one. By comparing square arrays and hexagonal arrays under the same ϕ, we find that there is no significant difference between two types of micropillar arrangement pattern. These results suggest that the arrangement pattern of micropillar does not have significant influence on the permeability of the arrays.

5 Conclusions

In this paper, we present a DPD model to simulate the flow through periodic arrays of circular micropillar and investigate the permeability of two types of micropillar arrays. Numerical results of permeability show a quantitative match with the theoretical solutions. First, we use the DPD model to simulate the flow around a single micropillar. The flow fields including the horizontal and vertical velocity fields, the number density field, and the streamline of the flow are analyzed. We also analyze the drag coefficient for the micropillar at different Reynolds numbers. The drag coefficient of circular cylinder is close to the experimental results for Re from 0 up to 100. These results indicate that the presented DPD method is useful and effective for the simulation of flow around a single micropillar. Later, we use the DPD model to sim¬ulate and investigate the permeability of two types of micropillar arrays. The predicted solid volume by the DPD simulation of two regularly distributed arrays is quite close to the actual counterparts. By comparing the two types of micropillar arrangement patterns, we find that there is obviously no difference in permeability of two types of micropillar arrays.