Advances in Computational Fluid Mechanics in Cellular Flow Manipulation: A Review

: Recently, remarkable developments have taken place, leading to signiﬁcant improvements in microﬂuidic methods to capture subtle biological e ﬀ ects down to single cells. As microﬂuidic devices are getting sophisticated, design optimization through experimentations is becoming more challenging. As a result, numerical simulations have contributed to this trend by o ﬀ ering a better understanding of cellular microenvironments hydrodynamics and optimizing the functionality of the current / emerging designs. The need for new marketable designs with advantageous hydrodynamics invokes easier access to e ﬃ cient as well as time-conservative numerical simulations to provide screening over cellular microenvironments, and to emulate physiological conditions with high accuracy. Therefore, an excerpt overview on how each numerical methodology and associated handling software works, and how they di ﬀ er in handling underlying hydrodynamic of lab-on-chip microﬂuidic is crucial. These numerical means rely on molecular and continuum levels of numerical simulations. The current review aims to serve as a guideline for researchers in this area by presenting a comprehensive characterization of various relevant simulation techniques.


Microfluidic Systems for Cellular Flow Manipulation
Most of the in vivo and in vitro studies have been carried out at a bulk-level. Although these techniques bring invaluable critical insights into various fields, including biology and medicine, they usually lack single-cell resolution features to recognize cell heterogeneities associated with diseases. For studying rare cell events and the intermediate extracellular signaling, through which cells communicate, analysis tools with the micro-level resolution are of significant importance. Microfluidics chips offer new approaches for cell assays and have also been used for studying of cell biology by providing platforms for manipulating, separating, sorting, filtering, trapping, and detecting tiny biological particles based on cellular heterogeneity. They can simulate small-scale fluid flow and chemical gradients and offer full manual control over the particles to study the desired details, e.g., for food market, clinical, pharmaceutical, and other applications [1][2][3][4][5][6][7]. Precise manipulations such as focusing, separation, and fractionation of cells is a vital capability of microfluidics which can be achieved by engineering hydrodynamics forces based on unique physical attributes of cells such as size [8,9], density [9,10], deformability [11][12][13], and morphology [14] using variety of methods such as crossflow filtration [15], electrode arrays [16], optical force switching [17] and other methods, some of and other methods, some of which can be found in detail in the review presented in [18]. In this spirit, hydrodynamic phenomena, which carry microenvironmental physical impacts on cells, are critical in almost all physiological functions and bodily systems. Generally, microfluidic devices can sort/isolate or partition cells by active approaches with external forces from acoustics [19,20], magnetics [21,22], and electrokinetics [23], as examples, as well as by passive approaches such as cellular immobilization [24], deterministic lateral displacement [25], and hydrodynamic filtration [26]. These devices have also been extensively used and proposed to study mechanical shear stress, and cell deformability and the recent advances of these applications can be found in [27][28][29][30]. This is extremely helpful since deformability of RBCs (red blood cells) can represent a potential sign of several sicknesses [28]. Moreover, microfluidic devices have been proposed to study the possible effects of microbubbles in microvessels [28] to represent the effects of these bubbles, which can form in the blood vessels occasionally and cause possible pathological events, such as preventing the food to reach to the cells in specific regions. Also, targeted drug delivery via nanomaterials, such as carbon nanotubes, was introduced recently and significantly increased the efficacy of the drugs [31][32][33][34]. As evaluation of these particles are highly dependent on their microenvironments, finding the optimum design point is facing some difficulties [35] and microfluidic systems can offer exceptional abilities to screen the nanoparticles to resolve these difficulties [1]. Zhu et al. [36] presented an in-dept review of nanoparticle delivery steps in cellular flows and summarized several microfluidic models specialized in nanoparticle evaluations related to such investigations. Particle sorting mechanisms are required to direct the drugs to the targeted cells, and numerical methods have proved themselves to be efficient in conducting such simulations [37,38]. Microfluidic devices, generally, implement any one of the two different approaches; active approaches and passive approaches. In the latter, as opposed to former which relies on fluorescent labels or beads, a variety of techniques based on inherent differences in cellular morphology among cell groups (e.g., size, shape, compressibility, and density) are employed to sort cells. In active manipulation approaches, based on specified field, the cell array is immersed in a space wherein energy distribution results in controlled forces that move cells along desired paths. On the other hand, in passive techniques, the two predominant opposite forces, known as wall lift force, which tends to repel cells from wall, and shear gradient lift force, which is due to the curvature of the fluid parcel velocities and tends to repel cells from center, determine cell equilibrium positions in the cross section of microchannel [39]. A summary of forces exerted on cells/particles is shown in Figure 1 for passive separation technique. It is noteworthy to mention that complex hybrid approaches have recently found their ways to microfluidic devices, and have been advantageous in certain studies [40]. Note that there is another force, Saffman force, which is channel center/wall-directed lift force experienced by particles lagging/leading the fluid and it is not shown in the figure to avoid confusions. Qualitatively, in most of the literature, that relied on passive techniques, and in a fair amount of the literature, focused on active approaches, the numerical simulation had been used to refine the Qualitatively, in most of the literature, that relied on passive techniques, and in a fair amount of the literature, focused on active approaches, the numerical simulation had been used to refine the understanding of physical phenomena, or/and optimize micro-channel geometrical aspects [41][42][43][44][45][46][47][48][49]. It is noteworthy that herein, we focus mostly on separating and trapping techniques on mammalian cells that are passive and suspended inflow. Hydrodynamics of swimming microorganisms such as bacteria, or adherent cell are not within this review's scope.

Modeling of Cell Separation/Partitioning
While utilizing simplified analytical analysis (such as Stokes flow) and trial and error-driven experimental approaches are prone to many unresolved issues, numerical simulation techniques are becoming a routine interrogation tool to avoid such problems. More importantly, numerical analyses are vital to unravel underlying physics, to test different hydrodynamic effects, and to optimize process conditions and proposed designs, which can be very expensive and require a long series of experiments or almost impossible experimentally with equipment at hand regardless of the scale of the experiments [50][51][52][53]. These approaches characterize microfluidic cell separation devices at two-particle and continuum levels. Formerly, numerical simulations were carried out to identify particle-level interactions using statistical mechanics over individual particle mechanics force field. In the latter, continuum-level simulations based on conservation laws are used to address microfluid system attributes. One can refer to lattice Boltzmann methods (LBM) as the main molecular level numerical simulation, and computational fluid dynamics (CFD) based on Navier-Stokes (NS) equations as the prominent continuum level simulation techniques utilized for cell separation applications ( Figure 2). In the following, fundamentals of these two numerical methods are briefly reviewed to convey the essential knowledge needed for the remainder of the work.

Modeling of Cell Separation/Partitioning
While utilizing simplified analytical analysis (such as Stokes flow) and trial and error-driven experimental approaches are prone to many unresolved issues, numerical simulation techniques are becoming a routine interrogation tool to avoid such problems. More importantly, numerical analyses are vital to unravel underlying physics, to test different hydrodynamic effects, and to optimize process conditions and proposed designs, which can be very expensive and require a long series of experiments or almost impossible experimentally with equipment at hand regardless of the scale of the experiments [50][51][52][53]. These approaches characterize microfluidic cell separation devices at twoparticle and continuum levels. Formerly, numerical simulations were carried out to identify particlelevel interactions using statistical mechanics over individual particle mechanics force field. In the latter, continuum-level simulations based on conservation laws are used to address microfluid system attributes. One can refer to lattice Boltzmann methods (LBM) as the main molecular level numerical simulation, and computational fluid dynamics (CFD) based on Navier-Stokes (NS) equations as the prominent continuum level simulation techniques utilized for cell separation applications ( Figure 2). In the following, fundamentals of these two numerical methods are briefly reviewed to convey the essential knowledge needed for the remainder of the work.

Lattice Boltzmann Method
When LBM first appeared on the front page of Washington Post on November 19, 1985, a dramatic enthusiasm was brought into computational methods, which was diminished later by revealed limitations, mostly because of convergence issues and computational expenses at the time, but these issues have been alleviated over the period since then. This CFD method is based on the kinetic theory, which correlates macroscopic properties of a medium to basic mechanical laws governing the behavior of its constituents as hard-sphere particles with elastic collision. State of each molecule is described by two variables, position(x(x,y,z)) and velocity (v(vx, vy,vz)) in space. Since it is challenging, if possible, to go after each molecule, instead, a velocity distribution function, f(x,v,t), which is the velocity number density of molecules at position x and speed v at a time t is defined; and then is used to obtain macroscopic flow properties such as pressure, temperature, density, etc. Then Boltzmann equation, which expresses a balance between particle/molecule transport (streaming) and collision must be analytically/numerically solved to find the mentioned velocity distribution function. This equation (in the absence of external forces) can be written as:

Lattice Boltzmann Method
When LBM first appeared on the front page of Washington Post on November 19, 1985, a dramatic enthusiasm was brought into computational methods, which was diminished later by revealed limitations, mostly because of convergence issues and computational expenses at the time, but these issues have been alleviated over the period since then. This CFD method is based on the kinetic theory, which correlates macroscopic properties of a medium to basic mechanical laws governing the behavior of its constituents as hard-sphere particles with elastic collision. State of each molecule is described by two variables, position(x(x,y,z)) and velocity (v(v x , v y ,v z )) in space. Since it is challenging, if possible, to go after each molecule, instead, a velocity distribution function, f(x,v,t), which is the velocity number density of molecules at position x and speed v at a time t is defined; and then is used to obtain macroscopic flow properties such as pressure, temperature, density, etc. Then Boltzmann equation, which expresses a balance between particle/molecule transport (streaming) and collision must be analytically/numerically solved to find the mentioned velocity distribution function. This equation (in the absence of external forces) can be written as: where x, v, and t represent position, velocity of the particles, and the associated time during the simulation. In this equation, C models the pairwise collision among particles. There exist, several models, to describe the function C: BGK, the classic and reliable method to solve problems with low Reynolds number (Re) (for more technical details, refer to [55]). MRT, enhances stability by presenting relaxing parameters to ensure stability [56]. A regularized model, renders better accuracy and stability compared to BGK by eliminating higher-order, non-hydrodynamic terms from the particle populations [57]; and Entropic model in which relaxation time is tailored locally to avoid an entropy decrease during the collision and ensuring unconditional numerical stability [58]. Another essential key input is the lattice types on which these models are operated. Several different lattices, both cubic and triangular, and with or without rest particles in the discrete distribution function exist. A typical way of clustering the various methods by lattice is the DnQm scheme. Here "Dn" stands for "n dimensions," while "Qm" stands for "m speeds." For example, D3Q15 is a 3-dimensional lattice Boltzmann model on a cubic grid, with rest particles present. Each node has a crystal shape and can deliver particles to 15 nodes: each of the 6 neighboring nodes that share a surface, the 8 neighboring nodes sharing a corner, and itself. It has been shown that the solution of the Boltzmann equation can represent the NS equations at low Mach number. Having this function, one can configure the 3D flow field of the system. Among many advantages, one can mention clear physical pictures with the detailed resolution, easy implementation of intricated boundary conditions with simple mechanical rules such as bounce-back and reflection and dynamic interfaces, and fully parallel algorithms for implementing high-performance computing. Nevertheless, as a natural-born dynamic scheme, the LBM is not a method of choice for steady-state analysis; and also is not well-versed for body-fitted coordinates and adaptive time-stepping [59]. As of now, there are several open-source software such as J-lattice-Boltzmann, SunlightLB, OpenLB, LIMBES, LBSim, LB2D, LB3D, Palabos, LBM-C, Taxila LBM, and Sailfish; and commercial ones like PowerFLow, XFlow, FlowKit, LBHydra and utraFluid taking advantage of this method. A schematic and step by step comparison between LBM and conventional CFD algorithms can be seen in Figure 3 below for a better understanding of the procedure in each numerical method.
adaptive time-stepping [59]. As of now, there are several open-source software such as J-lattice-Boltzmann, SunlightLB, OpenLB, LIMBES, LBSim, LB2D, LB3D, Palabos, LBM-C, Taxila LBM, and Sailfish; and commercial ones like PowerFLow, XFlow, FlowKit, LBHydra and utraFluid taking advantage of this method. A schematic and step by step comparison between LBM and conventional CFD algorithms can be seen in Figure 3 below for a better understanding of the procedure in each numerical method.  This method has been massively used in simulations of RBCs to study cell mechanical deformation behaviors under certain circumstances and/or suspension of these particles in various channel shapes representing the blood vessels to help the sorting application of microfluidic devices. A few examples of how LBM method performed in the real applications are presented in the following. Deformability of RBCs has been studied using LBM in simple channel conditions [60,61] and/or under more complex environments [62,63]. A few examples of results obtained from simulations using LBM method to simulate phenomena related to RBCs are presented in Figure 4. For a good review of such simulations in the past and more examples, one can refer to [64], where the researchers have presented results obtained from LBM simulations of RBCs-related phenomena.
A few examples of how LBM method performed in the real applications are presented in the following. Deformability of RBCs has been studied using LBM in simple channel conditions [60,61] and/or under more complex environments [62,63]. A few examples of results obtained from simulations using LBM method to simulate phenomena related to RBCs are presented in Figure 4. For a good review of such simulations in the past and more examples, one can refer to [64], where the researchers have presented results obtained from LBM simulations of RBCs-related phenomena.

Conventional CFD Methods
These methods are based on continuum assumption and include discretizing conservation equations coupled with boundary conditions. Widely accepted conventional CFD techniques include finite difference, finite volume, and finite element methods (FEM). There are other, yet infrequent methods, such as boundary element, spectral element, and other in-house developed software tools. There are two reviews covering these tools in the context of microfluid devices in general and, capturing rare cell systems. Erickson [69] provided an overview of traditional and emerging CFD techniques for integrated microfluidic devices covering thermal and chemical microscale flow and

Conventional CFD Methods
These methods are based on continuum assumption and include discretizing conservation equations coupled with boundary conditions. Widely accepted conventional CFD techniques include finite difference, finite volume, and finite element methods (FEM). There are other, yet infrequent methods, such as boundary element, spectral element, and other in-house developed software tools. There are two reviews covering these tools in the context of microfluid devices in general and, capturing rare cell systems. Erickson [69] provided an overview of traditional and emerging CFD techniques for integrated microfluidic devices covering thermal and chemical microscale flow and species transport simulation and cellular/particulate transport. On the other hand, the recent work of Jarvas and Guttman [70] centers on overviewing these tools on cell sorting and rare cell capturing, such as circulating tumor cells. Herein, a brief introduction of some of the main traditional CFD tools is given. These techniques can be clustered by the particular meshing approach to discretize the governing equations such as finite difference, finite volume, FEM. In finite difference method, the governed partial differential equations are approximated with Taylor series based on the values of neighboring nodes. This, then, renders a system of algebraic equations to be solved using a variety of well-developed explicit and implicit solvers. This method has the advantage of generally being the easiest to implement, though, it is practically limited to well-structured grids and thus, only advantageous in simple geometries (e.g., capillaries or channel cross-sections) (refer to Taflove and Hagness [71] for more detail). Writing codes mostly use this method in generic programming software such as C, C++ FORTRAN, etc.
In finite volume method, partial differential equations describing conservation laws are integrated and evaluated at the surfaces and finally discretized for very small, but finite-sized elements constituting the whole body. Unlike finite difference method, this method does not demand a structured grid and is thus, suitable for more complex and commonly encountered geometries (e.g., looping channels used in on-chip capillary electrophoresis). Since the conservation equations are to be applied on each irregular volume, the definition of derivatives and following Taylor series expansion to represent such is very challenging. It has been shown that this method is best suited for problems where the viscous terms are absent as opposed to the low Reynolds number flows encountered in microfluidics, in which these terms are dominant [71]. This method is incorporated in many commercialized software packages. For more detail on fundamentals, refer to [72].
In FEM, like finite volume method, the whole body is subsidized to smaller, simpler parts that are called finite elements. The simple equations that model these finite elements are then systematically recombined to form a more extensive system of equations that models the entire elements. Unlike the finite volume method, at this point, whereby the PDEs are to be discretized, in FEM variational methods from the calculus of variations are used to approximate a solution by minimizing an associated error function. FEM shares many of the same advantages of the finite volume method, most importantly, the ability to manage arbitrary grids and complex geometries.
From the microfluidicist's outlook, although it is mathematically expensive and laborious, in terms of implementing complex boundary conditions, with gradient-based boundary conditions that are commonly encountered in applied electrical fields, transport systems with surface phase reactions and thermal analysis involving convective heat transfer, this method is much easier to work with compared to the aforementioned ones [73]. Isoperimetric quadratic elements can also be employed to precisely conform to curved fluid-fluid interfaces such as those encountered in capillary flows [69]. Among others, the main disadvantage of FEM is the mathematical difficulties associated with handling highly irregularly shaped elements (e.g., large aspect ratios or highly skewed) or large ranges in element size within a single mesh. In lab-on-a-chip devices, the relevant length scales can range over seven orders of magnitude, from the double-layer thickness (nm) to channel length (cm) and thus, it is often difficult to avoid using such elements while maintaining a computationally tractable problem. Despite these challenges, dramatic developments in computational capabilities during recent years make this method pivotal for several well-known commercial software packages such as NASTRAN, ANSYS, ABAQUS, COMSOL, and FIDAP.
A few examples of how these methods performed in the real applications are presented in the following. CFD methods have mostly been applied to investigate the phenomena related to human arteries, such as blockage because of atherosclerosis, the application and performance of ventricular assist devices (VADs) and etc. Song et al. [74] implemented CFX to study the performance of a CF4 type VAD pump and effects of the pump on the flow path. Chesnutt et al. [75] used the discrete element method to study transport, activation, and adhesion of RBCs through thrombus formation around endothelium. Jung and Hassanein [76] used a three-phase CFD method to investigate RBCs reflections in disturbed flow regions of human arteries. Another interesting study was conducted by [77], where the authors studied the dynamics of RBCs in capillaries of finger nail-fold for various cases. A more recent study was done by [78], where the researchers used the CFD discrete element method (CFD-DEM) to conceptually design a microfluidic device capable of magnetic sorting of malaria-infected red blood cells. A few examples of results obtained from simulations using CFD methods in similar applications are presented in Figure 5. More details about applications of these methods can be found in [69,[79][80][81].
Appl. Sci. 2019, 9, x FOR PEER REVIEW 7 of 23 cases. A more recent study was done by [78], where the researchers used the CFD discrete element method (CFD-DEM) to conceptually design a microfluidic device capable of magnetic sorting of malaria-infected red blood cells. A few examples of results obtained from simulations using CFD methods in similar applications are presented in Figure 5. More details about applications of these methods can be found in [69,[79][80][81].

Scope of This Work
This review reports the use of computational fluid modeling techniques to develop and optimize new passive and active microfluid systems for analysis and control of cellular transport phenomena. Capabilities of these numerical methodologies along with their merits and cons will be extensively detailed so that multidisciplinary readers would be able to decide which method meets their requirements. Nevertheless, this article is written assuming readers have at least a rudimentary understanding of numerical simulation procedure. As the range of physical phenomena covered in this article is so extensive and the number of different techniques used to model each situation varies so greatly, a comprehensive review of the governing equations is beyond the scope of this article. Here, the focus is on the applications, the specific numerical tools, and the approaches used to model the relevant physics. Readers will be referred to relevant literature covering specifics and fundamental concepts.

Outlook
Lack of standard guideline for using CFD tools, either commercial or generic, was found in the literature, which initiated the present work. As microfluidic systems continue to mature and are applied to increasingly complex applications, we will continue to see their use in cancer and malaria studies and other diseases in general. Microfluidics is well suited to cell isolation, cell culture, threedimensional culture, cell/tissue perturbation, and analysis. Integration of these functions and increasing use of co-culture methods will undoubtedly play a significant role in cancer analysis in the coming years. From biochemical to mechanical assays, microfluidic chips are uniquely poised to make high impact discoveries, fueled less by innovations in device design and more by innovations in the questions asked with those devices. As mentioned before, an inexpensive way of assuring the

Scope of This Work
This review reports the use of computational fluid modeling techniques to develop and optimize new passive and active microfluid systems for analysis and control of cellular transport phenomena. Capabilities of these numerical methodologies along with their merits and cons will be extensively detailed so that multidisciplinary readers would be able to decide which method meets their requirements. Nevertheless, this article is written assuming readers have at least a rudimentary understanding of numerical simulation procedure. As the range of physical phenomena covered in this article is so extensive and the number of different techniques used to model each situation varies so greatly, a comprehensive review of the governing equations is beyond the scope of this article. Here, the focus is on the applications, the specific numerical tools, and the approaches used to model the relevant physics. Readers will be referred to relevant literature covering specifics and fundamental concepts.

Outlook
Lack of standard guideline for using CFD tools, either commercial or generic, was found in the literature, which initiated the present work. As microfluidic systems continue to mature and are applied to increasingly complex applications, we will continue to see their use in cancer and malaria studies and other diseases in general. Microfluidics is well suited to cell isolation, cell culture, three-dimensional culture, cell/tissue perturbation, and analysis. Integration of these functions and increasing use of co-culture methods will undoubtedly play a significant role in cancer analysis in the coming years. From biochemical to mechanical assays, microfluidic chips are uniquely poised to make high impact discoveries, fueled less by innovations in device design and more by innovations in the questions asked with those devices. As mentioned before, an inexpensive way of assuring the applications of these devices is implementing the numerical methods developed for such purposes. In the following, a brief summarization of these various methods is provided and categorized based on their characteristics. A brief illustration of these methods, based on the appropriate problem dimensions that they can be used for, is presented in Figure 6 below. This graph provides a simple benchmark for choosing the appropriate numerical model according to the problem-specific characteristics and dimensions. applications of these devices is implementing the numerical methods developed for such purposes. In the following, a brief summarization of these various methods is provided and categorized based on their characteristics. A brief illustration of these methods, based on the appropriate problem dimensions that they can be used for, is presented in Figure 6 below. This graph provides a simple benchmark for choosing the appropriate numerical model according to the problem-specific characteristics and dimensions. Figure 6. Various numerical models for simulation of appropriate particle problems based on the application problem dimensions. Adapted from [59].

LBM
As mentioned earlier, LBM is considered an efficient solver for the NS equations for microfluidic systems because of its particulate-based nature and parallel computing advantages. LBM is typically viewed as a second-order accurate numerical method in space and time. Since this method emerged quite recently, the use of it for cellular fractionation purposes is relatively limited. In simulating blood cells or deformable particles immersed in a fluid, though, a higher number of articles have been published. Because of similarity, the configuration of these problems and the numerical setup can be used for investigating microfluidic systems. The lattice Boltzmann method has been massively used to study the cellular flow and related topics in the recent literature and found extremely useful to deal with difficulties of such problems.
LBM has several advantages over other conventional CFD or particle-based methods and has attracted a lot of interests in computational physics. Because of its particle-based nature, it is also easy and flexible to model complex structure fluids, like DPD and SPH. As a mesoscopic method, it is more straightforward to incorporate microscopic interactions, like DPD. Besides, the primary advantage may be that it can be easily implemented in a massive parallel computing environment because of its local dynamics. Finally, it deals with the fluid-RBC interaction by coupling to the immersed boundary method, instead of modeling the RBC membrane as physically connected particles, which is entirely different. Comparing with the conventional methods for multiphase flows, the LBM does not track interfaces, while sharp interfaces can be maintained automatically.

LBM
As mentioned earlier, LBM is considered an efficient solver for the NS equations for microfluidic systems because of its particulate-based nature and parallel computing advantages. LBM is typically viewed as a second-order accurate numerical method in space and time. Since this method emerged quite recently, the use of it for cellular fractionation purposes is relatively limited. In simulating blood cells or deformable particles immersed in a fluid, though, a higher number of articles have been published. Because of similarity, the configuration of these problems and the numerical setup can be used for investigating microfluidic systems. The lattice Boltzmann method has been massively used to study the cellular flow and related topics in the recent literature and found extremely useful to deal with difficulties of such problems.
LBM has several advantages over other conventional CFD or particle-based methods and has attracted a lot of interests in computational physics. Because of its particle-based nature, it is also easy and flexible to model complex structure fluids, like DPD and SPH. As a mesoscopic method, it is more straightforward to incorporate microscopic interactions, like DPD. Besides, the primary advantage may be that it can be easily implemented in a massive parallel computing environment because of its local dynamics. Finally, it deals with the fluid-RBC interaction by coupling to the immersed boundary method, instead of modeling the RBC membrane as physically connected particles, which is entirely different. Comparing with the conventional methods for multiphase flows, the LBM does not track interfaces, while sharp interfaces can be maintained automatically.
For instance, it is easy to incorporate mesoscale physics, such as interfacial breakup or coalescence. Moreover, the computational cost for simulating realistic fluid flows is reasonable compared with particle-based methods (e.g., molecular dynamics). Also, the physics associated with molecular-level interactions can be incorporated more quickly, since the Boltzmann equation is kinetic-based. Hence, the lattice Boltzmann model might be fruitfully applied to micro-scale fluid flow problems to some extents. Using Lagrange multiplier as a coupling term for fluid-solid interactions in LBM corresponds to calculating the displacements of the solid particles implicitly, which is proven to be more precise and accurate compared to immersed boundary method, which makes it a favorable method in the modeling of RBC deformation modeling. Moreover, this method offers other advantages over other conventional CFD methods for simulation of particle suspension such as the possibility of deriving different simplified local particle-fluid interaction rules because of its direct connection to molecular mechanics [84].
The first complication, and often most important numerically, is the range of relevant length scales, which can vary as much as seven orders of magnitude (for instance, from the double layer thickness, nm, to channel lengths and substrate dimensions, cm). For the second limiting factor, one needs to modify the normal lattice Boltzmann method to account for the particles existing in the fluid and their interaction as boundary conditions in the method. Researchers have modified the LBM method and implemented in studying different phenomena related to cellular and blood flow study to overcome some of these limitations. [85] implemented a red blood cell model (RBC) in an LBM solver to investigate and model red blood cells and fluid motion. For this purpose, they employed an LBM solver to an immersed boundary scheme to study the fluid motion. In the frame of this modified method, otherwise called IB-LBM, Lagrangian force, F(s,t), acting on the fluid by each finite segment of the domain boundary is calculated based on the boundary configuration near that and then is distributed to the ambient fluid in the vicinity of the boundary. Thus, treating it as a body force (refer to [86] for more detail) and then the fluid motion can be solved using the lattice Boltzmann equation and position(x(x,y,z)) and velocity (v(v x , v y ,v z )) in space can be updated. Iterating over the time and updating position and velocity functions after each iteration, the simulation continues until the convergence conditions are met. When the particles' density distributions are found, the macroscale quantities can be calculated using different averaging methods. [87] coupled lattice Boltzmann method with the RBC model and modeled the fluid domain with a simple external force. In this method, the RBC network is used to model the red blood cells deformation and lattice Boltzmann method is employed to track the particles and compute the averaged quantities over time. This method is used to solve a domain including fluid and solid subdomains having the Lagrange multiplier served as the external force representing a coupling term for solid and fluid subdomains. Another interesting modification of this method was conducted by [88], where the authors coupled LBM and FEM methods using immersed boundary method to investigate the behavior of RBCs in rectangular channels under different conditions to represent the cells' behavior in microfluidic devices. In this study, they studied RBCs membrane mechanics by using FEM using Cauchy stress tensor and defining virtual displacement and stress tensor, while solving for fluid mechanics using LBM.

DPD
Dissipative particle dynamics method (also called DPD) is a Lagrangian, and coarse-grained mesoscale particle-based hydrodynamic method, in which the domain is represented by a distribution of discrete particles with separated physical properties and each particle consists of a collection of atoms and molecules, that was initially introduced to deal with complex Newtonian fluids and soft matters and is ideal for cell biophysics problems, such as blood flow behavior and blood flow particle interactions [89][90][91][92]. Having a particle-based nature along with being able to reach higher time scales than methods such as MD, make it a relatively valuable and computationally efficient method for microfluidics in terms of accuracy in capturing the details and capability to reach higher time scales [90]. As discussed before, continuum methods have been proved to be applicable for many of the small-scale simulations, however, for some complex problems, there is a need for more microscopic details out of their grasps. There are many cases, where given a large number of particles such as proteins and cells in the computational domain, high computer resources are demanding. In this situation, applications of continuum-based models, where subtle microscopic delicacies and their interactions with the mainstream are ignored, might pose further precision issues. Mesoscale methods can, thus, be useful for such problems. DPD is clustered under these methods and has proved itself an appropriate method for a variety of problems with significant hydrodynamics details and/or thermal fluctuations. DPM can be easily implemented to the code by just changing the equation used to model the conservative forces among the particles [90]. This method, as discussed, is capable of reaching higher time scales than some other particle-based method such as MD. However, the accuracy of the method depends highly on the time step chosen for simulation, and thus, this can be a tricky part specific to the problem in hand [90].
Original DPD method was, however, too simplistic, for instance in modeling the friction forces, and thus, not appropriate to capture the real physics in detail in many complex problems consisting of different types of particles. To maintain the accuracy of DPD method, one needs to input real values for some physical parameters, such as viscosity, which is case dependent, i.e., one would need to extract real data from experiments conducted on a comparable problem with the same details before using the method to numerically investigate the same problem [90,91]. Lack of freedom of using different equations of states during the simulation makes the DPD method susceptible to having unrealistic thermodynamic behavior for complex systems, including complex simulations of blood flow and sustaining temperature gradients. In particular, for higher gradients in temperatures, because of lack of energy conservation in the method, unreliable results might be anticipated [90]. Furthermore, this method has yet deficiencies to face to model unbounded particle interactions [90].
There have been many improvements proposed in the literature since the arrival of the original method to make it capable of dealing with more complex fluids and systems and adopting with newly developed models to capture more accurate details in fluid flow simulations. The most important modified versions of this method are many-body DPD (MDPD), energy-conserving DPD (EDPD), fluid particle method (FPM) and finally, smoothed dissipative particle dynamics method (SDPD) and each tackled some of DPD drawbacks. For instance, MDPD introduced an independent equation for temperature distribution and a more general equation of state to ensure having more realistic temperature distribution, EDPD added an extra internal energy value for every particle for non-isothermal cases, and FPM tackled the problem of having too simplistic friction forces in the original method. Finally, SDPD, which is derived from SPH and is similar to MD in terms of moving particles interactions simulation, introduced Lagrangian discretization of NS equations and an exact discretized fluctuation theorem and adopted a thermostat similar to the one in DPD, improving most of the original method's drawbacks [90]. This method has all the advantages of the other three methods and none of their disadvantages, which made it the most common DPD method. SDPD has several advantages over the conventional DPD method, such as allowing direct input of fluid transport coefficients and selecting desired arbitrary equation of state; resulting in full control over fluid compressibility [91]. This method, however, has all the drawbacks of the SPH method because of the closely related particle-based nature of it compared to SPH [90]. Thus, one would still need to input high accuracy real physical values, such as viscosity, to the method to reduce the inaccuracies. Moreover, the pressure field of the problem needs to be calculated separately by using the equation of state via solving the Poisson equation of pressure, given its inherent continuum-originated build-in models [93]. Another drawback relates to the lack of a direct relationship or link between the viscosity of the fluid and the velocity profile, which was alleviated in SDPD [94].
These methods have been extensively used for simulation of particle sorting systems, for instance, in deterministic lateral displacement devices, modeling blood flow problems including RBCs related simulations, blood-related diseases, and disorders, particles like plasma and proteins, bacteriology and parasitology over the recent years [89,[91][92][93][94]. SDPD method has proven to be especially suitable for simulation of blood flow and by being coupled with the immersed boundary method. The fluid-RBCs interactions, which are an essential area of focus in cell sorting and cell separation problems, can also be accurately modeled instead of being modeled as physically connected particles, which makes it an accurate and computationally efficient method for blood flow problems [93]. This method can also be applied to understand the behavior of healthy and unhealthy RBCs in real physics as well.
For years, lacking appropriate boundary conditions and high-performance computing hindered implementing particle-based methods in simulation of blood flow in microvascular networks and similar complex problems. Likewise, the SDPD method suffered the same hindrances and could not be used for many of the investigations, despite all other improvements since the time the original DPD method was introduced [95]. Some issues to be tackled in terms of the computational speed of this method are; for instance, irregular patterns of computation resulting in hindering its parallel computing capability and the maximum speed limit and having an order of magnitude higher interaction instructions than in MD for each time step. Building a neighborhood list in each time step for each of the particles and finally, implementation of an extra thermostat for some of the forces and interactions [95]. New computers and algorithms for higher performance computing were presented recently and are shown to solve these two problems in the future. Recent advances in numerical computing capabilities and new parallel algorithms for these particle-based methods made the implementation of numerical simulations easier and capable of reaching higher time scales, which is a crucial factor for simulation to replace the expensive experimental tests on devices.
Aside from advances in computational technologies, efforts have been devoted to increasing the method's efficiency and speed and thus, reaching higher computation times and computational domains with a higher number of entities. For instance, one can mention the development of new open flow boundary conditions, which are clustered in the periodic boundary condition category. These boundary conditions still have some drawbacks for cases with so many inlets and outlets and also with time-dependent non-period flows, such as the system of the arterial network in human cardiovascular systems [95]. A few examples of results obtained from related simulations using DPD-originated methods are presented in Figure 7.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 11 of 23 connected particles, which makes it an accurate and computationally efficient method for blood flow problems [93]. This method can also be applied to understand the behavior of healthy and unhealthy RBCs in real physics as well. For years, lacking appropriate boundary conditions and high-performance computing hindered implementing particle-based methods in simulation of blood flow in microvascular networks and similar complex problems. Likewise, the SDPD method suffered the same hindrances and could not be used for many of the investigations, despite all other improvements since the time the original DPD method was introduced [95]. Some issues to be tackled in terms of the computational speed of this method are; for instance, irregular patterns of computation resulting in hindering its parallel computing capability and the maximum speed limit and having an order of magnitude higher interaction instructions than in MD for each time step. Building a neighborhood list in each time step for each of the particles and finally, implementation of an extra thermostat for some of the forces and interactions [95]. New computers and algorithms for higher performance computing were presented recently and are shown to solve these two problems in the future. Recent advances in numerical computing capabilities and new parallel algorithms for these particle-based methods made the implementation of numerical simulations easier and capable of reaching higher time scales, which is a crucial factor for simulation to replace the expensive experimental tests on devices.
Aside from advances in computational technologies, efforts have been devoted to increasing the method's efficiency and speed and thus, reaching higher computation times and computational domains with a higher number of entities. For instance, one can mention the development of new open flow boundary conditions, which are clustered in the periodic boundary condition category. These boundary conditions still have some drawbacks for cases with so many inlets and outlets and also with time-dependent non-period flows, such as the system of the arterial network in human cardiovascular systems [95]. A few examples of results obtained from related simulations using DPDoriginated methods are presented in Figure 7.   [96]. (b) DPD simulation of deformation of a single healthy RBC over going through a tube flow [94]. (c) Deformation of a single RBC going through a micro channel [96].
As this method has been improved over recent years, some commercial software packages have been developed using mainly FORTRAN and C++ to simplify the process and to create a benchmark for future improvements. Table 1 summarizes some of the most common commercial tools developed based on these models [97]. Note that some of these tools are capable of using different or multiple mesoscale algorithms to model a system of particles. Table 1. Some software packages developed based on the original DPD or the associated improved versions [97].

Tools
Multiple Algorithms These software tools are mainly developed using C++ and FORTRAN languages and have easy interfaces to be used by a typical user. However, one would need basic knowledge of the algorithm they are created based upon as well as a basic understanding of available libraries in Windows and Linux operating systems depending on the case under investigation. In the following, Figure 8 illustrates a schematic of the DPD algorithm, which imparts a better understanding of necessary steps for developing codes using this method.

Mesocite
Appl. Sci. 2019, 9, x FOR PEER REVIEW 12 of 23 As this method has been improved over recent years, some commercial software packages have been developed using mainly FORTRAN and C++ to simplify the process and to create a benchmark for future improvements. Table 1 summarizes some of the most common commercial tools developed based on these models [97]. Note that some of these tools are capable of using different or multiple mesoscale algorithms to model a system of particles. These software tools are mainly developed using C++ and FORTRAN languages and have easy interfaces to be used by a typical user. However, one would need basic knowledge of the algorithm they are created based upon as well as a basic understanding of available libraries in Windows and Linux operating systems depending on the case under investigation. In the following, Figure 8 illustrates a schematic of the DPD algorithm, which imparts a better understanding of necessary steps for developing codes using this method.

Other Particle-Based Methods
There are other particle methods to consider to investigate such problems numerically. Two of the most common ones are direct simulation Monte Carlo and molecular dynamics methods. These two methods are particle-dynamics based methods and are highly accurate, although they are very computationally expensive and are particularly applicable for problems with very small dimensions. These methods can, nevertheless, be implemented for all problems in this area. They both are governed by the Boltzmann equation provided below: The original MD method was introduced in the 1950s based on some scattered physical theories [98] and has been modified multiple times ever since. Many of the researchers implemented this method to study the possibility of different nanomaterials to be used in cancer treatment [32,[99][100][101]. In this method, a number of particles are randomly distributed in the computational domain and a velocity distribution formula is used to assign each of them with a specific velocity, which depends on the problem being investigated. Many different velocity distribution functions have been introduced, such as Maxwell distribution function [102], with specific characteristics and are available in the literature. From this point, the particles are assumed to be moving with their velocity at the beginning of each time step being affected and being subject to boundary conditions and the interaction laws among them. These interactions are modeled by pre-defined potential functions derived from the second law of Newton and form the pressure field and thus, the force field, which can be monitored in detail throughout the simulation [103]. Many different possible functions have been introduced according to the needs for the case being investigated and from them, one can mention Lennard-Jones (LJ) potential function, which is the most common function and mathematically simple to implement. Efforts have been devoted to improve the inherent characteristics of this method and make it more computationally efficient to implement this method to bigger length-scale problems as well. For such developments, one can mention the improvements and novel methods in computing the force and pressure fields, which will be addressed in the following. [104] introduced a new force field for CHARMM36m with improved accuracy in generating polypeptide backbone ensembles for complex and disordered proteins. A new indirect method to calculate the forces in MD methods was introduced by [105], which is similar to the method used in macro-scale simulation for calculating the average force exerted by the flow on any external bodies.
This method has offered several advantages over other conventional CFD or particle-based methods. This method is grid-free, meaning that to model a problem, equations need only to be satisfied with respect to the geometry of the domain control volume, which makes it easier to implement compared to CFD methods. This method, however needs a new grid geometry, called bins, to account for cut-off radius and to average the quantities in the domain to calculate the macroscopic values for each quantity. This advantage, its inherent high accuracy with implementing periodic boundary conditions allow simulation of flow in complex geometries, regardless of the number of particles. Moreover, equations for this method are independent of the number of particles. Thus, one can assume smaller control volumes in a big domain to speed up the process and simplify the mathematical procedure. A simplified schematic of the algorithm for the MD method is presented in Figure 9 below. DSMC is another particle method to be used to investigate cellular flow problems. This method was introduced by [106] and is similar to MD in tracking particle properties and averaging their quantities to compute the macroscopic properties. Since the algorithm is a lot similar to that of MD, it is not schematized separately in this work. For a more detailed summary of the recent improvement, one can refer to [107]. The most important advantage of this method over MD is higher efficiency in terms of computational requirements, which makes it the preferred method in many cases because of computational tools limitation. The difference between this method and MD is in the procedure treating the interaction between the particle throughout the computational domain and this, in turn, impacts the necessary time step to ensure that the physics can be captured accurately. In MD, this time step is determined by comparing the particles' trajectories, while in DSMC, a simplified statistical equation is used to compute this time based on some assumptions. This statistical equation makes the mathematical formulation of this method simpler to understand and implement in numerical codes, since determining the time step is independent of the geometry. However, these simplifications sometimes lead to less accuracy or not trustable results in more sensitive cases, especially in cases with larger dimensions. A practical comparison of these two methods is presented in [108]. Researchers have modified this method over time to improve its accuracy over time and presented modified versions of this method to be implemented for cellular flow problems [109].
Smoothed particle hydrodynamics (SPH) is another method that has been used to investigate immersed particles in fluids over the years. This method is considered macroscopic, but particlebased and was proposed by [110]. The original method was too simplistic, but, has been improved over time and found its way for simulations of higher scale problems [111,112]. This method is one of the oldest mesh-free methods and can model complex problems or geometries, which is one of the advantages of this method over CFD methods. The nature of this method is similar to that of DPD, in terms of the mathematical formulation and can also be compared to MD, where it is derived from second law of Newton. However, this method differs from DPD or SDPD methods since it starts with DSMC is another particle method to be used to investigate cellular flow problems. This method was introduced by [106] and is similar to MD in tracking particle properties and averaging their quantities to compute the macroscopic properties. Since the algorithm is a lot similar to that of MD, it is not schematized separately in this work. For a more detailed summary of the recent improvement, one can refer to [107]. The most important advantage of this method over MD is higher efficiency in terms of computational requirements, which makes it the preferred method in many cases because of computational tools limitation. The difference between this method and MD is in the procedure treating the interaction between the particle throughout the computational domain and this, in turn, impacts the necessary time step to ensure that the physics can be captured accurately. In MD, this time step is determined by comparing the particles' trajectories, while in DSMC, a simplified statistical equation is used to compute this time based on some assumptions. This statistical equation makes the mathematical formulation of this method simpler to understand and implement in numerical codes, since determining the time step is independent of the geometry. However, these simplifications sometimes lead to less accuracy or not trustable results in more sensitive cases, especially in cases with larger dimensions. A practical comparison of these two methods is presented in [108]. Researchers have modified this method over time to improve its accuracy over time and presented modified versions of this method to be implemented for cellular flow problems [109].
Smoothed particle hydrodynamics (SPH) is another method that has been used to investigate immersed particles in fluids over the years. This method is considered macroscopic, but particle-based and was proposed by [110]. The original method was too simplistic, but, has been improved over time and found its way for simulations of higher scale problems [111,112]. This method is one of the oldest mesh-free methods and can model complex problems or geometries, which is one of the advantages of this method over CFD methods. The nature of this method is similar to that of DPD, in terms of the mathematical formulation and can also be compared to MD, where it is derived from second law of Newton. However, this method differs from DPD or SDPD methods since it starts with NS equations and discretizes the computational domain into finite points that are treated as groups of particles using Lagrangian discretization. A good review of SPH method is presented in [113]. In addition to being a mesh-free method, unlike other particle-based methods, this method works directly with macroscopic quantities because of NS equation being the basis for the computational steps in SPH. SPH, however, is not capable of addressing the mesoscopic problems (as opposed to SDPD), where there are large thermal fluctuations introduced because of the lack of a fluctuation-dissipation equation/function. Given the premise of this method for Newtonian fluids, it is facing some issues when dealing with non-Newtonian flows. Because of its inherent drawbacks regarding mesoscale problems and mesoscale nature of RBC related problems, the application of this method has been very limited in simulation of RBCs problems. One may refer to the study conducted in [112], where the researchers studied the deformability of Malaria-infected RBCs using in house 3D code developed based on SPH method. A review comparing the applications of this method with those of MD is presented in [114].
There are other mesh-free methods available, which can especially be used for studying deformability and mechanical properties of RBCs undergoing effects of various microfluidic devices or some diseases, such as discrete particle (DP), element-free Galerkin method (EFG), moving particle semi-implicit method (MPS), moving particle finite element method (MPFEM), immersed particle method (IPM), boundary element method (BEM), and multi-particle collision dynamics method (MPCD). However, these methods have not attracted enough attention in this field, and not many studies using them can be found. Most of these methods are mesh-free, which makes them merely implementable for complex fluid-structure problems as mentioned before, and are similar, in nature, to one of the particle-based techniques discussed previously.
BEM is another mesh-free method, in which computational discretized points are only defined on the physical and inlet-/outlet boundaries of the problem, and the analytical solutions for velocity fields in the problem results from integration formulations on the boundary and the quantities inside the boundaries can be calculated using somewhat complex integration formulations [115] depending on the geometry and condition of problem being investigated. An advantage of this method is that it can be easily used for problems with complex geometries since there is only a computational grid required for the boundaries. By using advanced integration formula, it is proven to be the most useful and accurate for infinite length problems compared to other similar methods. However, this method is determined to be the most reliable for stokes flow regimes, and it is not applicable for some cellular flow problems with a length-scale exceeding 100 micrometers and other methods are needed to be chosen instead [116]. MPS is very similar to SPH method in its nature, where the governing equations are derived from NS equations and the physical properties are calculated using a similar function as the kernel function in SPH method [117] and using weighted averaging functions based on the distance between two particles. The mathematical implementation of MPS is somewhat difficult if trying to reach a higher accuracy, which can be considered a drawback by some researchers, while the accuracy is not noticeably different. DP is another type of these methods, which is a Meso-scale method similar to original DPD or SDPD, where different forces of different nature are defined and act separately on the particles, i.e., conservative, dissipative and random forces. The difference between this method and SDPD, however, is the existence of new function to control and take into account the rotation of each particle throughout the simulation.
Some of these methods have had minimal applications in the literature, and thus, an in-detail explanation of them is out of the scope of this work. A short review on some of these methods, however, can be found in [64,114,118]. A few examples of the results obtained from related simulations using these particle methods are presented in Figure 10. Note that, the application of simulations conducted based on the original forms of these methods is not common in RBC-related investigation, however, the hybrid versions of these methods with DPD or LBM methods can be seen in the literature as it was mentioned earlier in this work.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 16 of 23 Figure 10. A few applications of moving-particle semi-implicit (MPS) and smoothed particle hydrodynamics (SPH) methods in simulation of RBCs. (a) SPH particle discretization of a single normal shaped RBC [112]. (b) Snapshot of flexible RBCs suspension in in blood flow using MPS method. [119].
As this method improved over recent years, some commercial software packages have been developed based on these methods to simplify the process and increase the efficiency of using these methods in investigating desired particle flow related problems. Table 2 summarizes some of the most common tools developed based on these models. These software tools are mainly developed using C++ and FORTRAN languages and have easy interfaces to be used by a typical user. However, one would need basic knowledge of the algorithm they are created based upon as well as a basic understanding of available libraries in Windows and Linux operating systems. Figure 10. A few applications of moving-particle semi-implicit (MPS) and smoothed particle hydrodynamics (SPH) methods in simulation of RBCs. (a) SPH particle discretization of a single normal shaped RBC [112]. (b) Snapshot of flexible RBCs suspension in in blood flow using MPS method. [119].
As this method improved over recent years, some commercial software packages have been developed based on these methods to simplify the process and increase the efficiency of using these methods in investigating desired particle flow related problems. Table 2 summarizes some of the most common tools developed based on these models. These software tools are mainly developed using C++ and FORTRAN languages and have easy interfaces to be used by a typical user. However, one would need basic knowledge of the algorithm they are created based upon as well as a basic understanding of available libraries in Windows and Linux operating systems.

Summary and Conclusions
The technologies employed in microfluidic devices have become much more sophisticated and complex over the past decade because of the significant developments of these devices in understanding the phenomena in human body cells. The pace of these improvements in microfluidic devices has accelerated over the past few years. Experimental methods to optimize these devices and understand related phenomena have lost their popularity because of the high costs associated with them and challenges on the way and numerical simulation methods have become the subject of interest to understand cellular-related phenomena better to help to optimize these devices. These methods, however, suffer from high computational costs. Thus, researchers have devoted their efforts to improve these methods, and new modifications have been introduced to help alleviate this problem and make them more applicable to study real-world problems.
The current paper gave a comprehensive and up-to-date review of various numerical methods available to study the cellular flow problems, with the primary focus on three highly used methods in body cells simulation: CFD, LBM, and DPD, but other particle-based methods and summarized and compared with these three methods as well. Also, schematics of algorithms were presented to be a guideline and starting point of understanding the procedure of implementing each technique for researchers. Also, this review discusses the practice areas for each method to be used depending on the characteristics of the problem being investigated and strengths of each method. The methods discussed in this work are summarized in Table 3 along with a short description of their applications and their associated dimensional nature. To emphasize the importance of each method for applicable problems, a summary comparison is provided in the following to conclude the benefits and drawbacks of each technique presented in this work. Table 3. Summary of main numerical methods presented in this work and applicability of them in similar simulations.

Method
Applications in the Literature

CFD
Mostly used for simulations of human arteries and to study VAD applications SPH Not frequently used for cellular flow simulations, but some RBC related simulation can be found recently

LBM
Frequently used for simulations of RBCs and related fluid-solid interactions as combined with IBM DPD Frequently used to study RBCs characteristics, blood flows, and behaviors of infected RBCs

Monte Carlo
Has been used for simulation of blood vessels, blood flow, and related problems, but lost its popularity over recent years MD Has been used for simulation of flow with complex structures, such as blood flow dynamics, but its application is limited because of high computational load • CFD methods have proven to be applicable to investigate the phenomena related to human arteries. However, to their more significant scale nature, they cannot be implemented to study smaller-scale problems related to body cells.

•
Many commercial software packages are available to study problems with high accuracy. • LBM is of highly parallel nature and can easily be implemented in massively parallel computing systems because of its local dynamics. LBM also offers excellent capabilities to be combined with other methods such as IBM to model RBC-flow interactions and deformation of RBCs. This is entirely different from modeling of the RBC membrane in physically connected particles like in SPH and DPD methods. • LBM and DPD, offer capabilities to model thermal fluctuation in the simulations, which is necessary to study phenomena such as RBCs aggregation. • LBM and DPD, as mesoscopic methods, can be used to study cellular flow problems with reasonable computational costs compared to microscale models. • SPH, unlike LBM and DPD, does not take into account the thermal fluctuations because of the lack of a fluctuation-dissipation function in its nature. However, this method originates from NS equation, and hence, one can obtain the physical properties of the flow directly using this method unlike the microscopic techniques such as MD, in which one needs to average the microscopic properties to get the physical macroscopic quantities.

•
Mesh-free methods such as SPH and other particle-based methods offer this capability to quickly model complex geometries without requiring complex computational grids as in CFD methods.
Finally, commercial packages available using the techniques discussed in this review were listed to cover the valuable and already available tools to researchers. The current review aims to serve as a guideline for scientists and researchers by gathering and summarizing techniques and their characteristics in one place. All in all, considering the advances in computational tools are becoming available for researchers in this field, particle-based methods and especially those of mesoscale nature are attracting more attention, and they need to be studied extensively and improved to be viable tools in the future to study large-scale cellular flow problems.
Author Contributions: S.P., N.K.A., and M.A. conceived this paper. S.P. and M.A. drafted the original manuscript. M.A., S.P., and N.K.A. critically revised it for relevant intellectual contents and contributions to the field. M.A, S.P., and N.K.A approved the final version of the paper.
Funding: This research received no external funding.