Numerical issues in design of passive micro-mixers and its performance analysis using CFD

Advection dominant fluid flows and transport processes in microfluidic mixing devices pose noteworthy challenges in design related simulations as regards numerical errors. In the present work, an unsteady flow model with decay constant was utilized for testing passive micro-mixer simulations for numerical error. The simulations were performed using OpenFOAM (Open-source Field Operation And Manipulation), a finite volume based open source software. Three types of hydrodynamic boundary conditions representing surface engineered device walls (hydro-active patches) were investigated as design test cases. Our study focuses upon simulated effect of these conditions on concentration profiles from the point of view of numerical error. We show that the model can be analytically treated in part to establish a particular feature in the simulated concentration profiles is a numerical artefact. Further, it is shown that this error is rectified either by changing the method of solving discretized equation or by coarsening the mesh. This process is also sensitive to the type of hydrodynamic boundary conditions imposed. A generalized grid convergence test we carried out attests to the reliability of these micro-fluidics simulation results, despite the fact that they display numerical artefacts. The paper addresses this anomaly to a strong interplay between difference equation solving methods, mesh parameters and boundary conditions in the partial differential equation (PDE) model for micro-flows in passive micro-mixers. Our investigation thus suggests that grid independence tests, instead of being based only upon cell variables and mesh parameters, should also be informed by PDE boundary conditions and difference equation solving methods. Keywords— Micro-fluidics , micro-mixer, PDE, CFD, OpenFOAM

INTRODUCTION Several experimental and numerical modeling studies have been performed for the design of micro-fluidic devices in which various geometric designs are suggested to enhance the mixing efficiency in these devices. Due to wide applications of micro-fluidic devices in varied disciplines these studies are important in the development of Lab-onchip technology. Micro-mixers are widely used in chemical, biological, medical, and environmental applications [1], [2], [3], [4]. Micro-mixers are usually classified as active and passive devices depending on the presence of external disturbance effects. External energy is employed to enhance mixing efficiency in active micro-mixers while passive do not require such external energy [2], [3], [4], [5], [6]. Passive micro-mixers have been widely used and explored by researchers in which improved flow, transport and mixing behavior of these devices are studied [7], [8]. Micro-scale mixing is more challenging process due to laminar flow which is performed at very low molecular diffusivities for many chemical and biological solutions [9], [10]. Time scale reactions of the micro-flows in micro-mixers are important to understand diffusion process [11], [12]. Computational Fluid Dynamics (CFD) simulations are helpful to investigate fluid flow as well as solute transport at micro scales. The numerical simulations are supported mostly on molecular models or on continuum models, depending on the type of application [13]. The continuum model being based on Navier-Stokes equations are used to describe fluid flow, when the micro-channels are in the range of the micrometers [14], [15], [16], [17], [18], [19]. Flow fields obtained from simulations are used in certain objective functions which are further optimized leading to design cases of micro-mixing © 2019, IJSRPAS All Rights Reserved 11 devices [20]. Finite volume method (FVM) and finite element method (FEM) are two numerical techniques commonly used in such simulations. It is known that FEM suffer from stability as well as numerical diffusion errors however FVM do not produce high false diffusion effects [21], [22]. FVM based software is the best choice for such type of simulations. It is crucial, in design and analysis of micromixers, that the level of accuracy of a given simulation of realistic flows to be estimated and controlled [19]. Particularly reliable and unsusceptible estimates of mixing performance of micro-mixers would critically involve error, estimation and control. To achieve this, it is important to know the influence of discretization and iterative convergence on numerical uncertainty in simulated concentration profiles [23]. Other influences can also be the algorithms for solving difference equations and nature of partial differential equations (PDE) boundary conditions. The degree of numerical errors in many solutions are related with the discretization schemes that are used in the solution of the governing equations (i.e., the magnitude of truncation errors) as well as other numerical and physical parameters used, such as grid size, grid type, and diffusivity [24]. Liu applied several tests using different numerical schemes under various grid size, flow, and transport conditions for a threedimensional micro-mixer [25]. Godunov suggested preferring higher order schemes due to their lower numerical diffusion effects in his study [26]. One example of such problems is an artificial viscosity. The effects of difference equation solving methods (not limited to micro-flows) are often presented in terms of error estimators. For example, an error equation method has been proposed for such estimation in laminar flows [27], [28]. Such issues of accuracy of error estimators are checked by comparing the results with exact analytical solution that are available (e.g. creeping stagnation flow, planar jet) [29]. Thus algorithms for solving difference equations have been addressed as causes of numerical error but the specific role of boundary conditions in conjunction with them needs attention. Numerical simulations of pressure drop for three dimensional rectangular microchannels have been presented where they proposed an equation for Hagenbach's factor of rectangular channel based upon these simulations [30]. Wall properties are important in such type laminar flow in micro-channels. It is noted that flow fields like velocity and quantities like friction factor exhibit heightened sensitivity to mesh refinement near walls [30]. An assumption is made that grid convergence studies using generalized criteria like Grid Convergence Index (designed by Roache) must be adequate to address such issues. These numerical errors are currently addressed using grid convergence tests suggested by Roache [31], [32], [33]. This criterion is uniformly applicable across various fluid flow simulations to be used for engineering design. Most of the authors addressing numerical simulation errors in fluid dynamics apply these criteria though some report specific spurious numerical artefacts without commenting on the suggested test [34].
In particular passive micro-mixers design would require specific device geometry based inputs regarding influence of numerical error in supporting simulations. In view of this, we have developed a model for unsteady flow with decaying chemical species and investigated it for three different design cases of a wall surface engineered passive micro-mixer. Related micro-flow simulations were carried out focusing on the issue of numerical error and were implemented using the open source software OpenFOAM. OpenFOAM software has been used for the purpose of benchmarking in studies proposing novel algorithms for solving difference equations related to Navier-Stoke's fluid model. It also emphasizes extensibility and flexibility so it is aimed first and foremost as a research tool [35], [36]. The particular design cases bore different boundary conditions of no-slip, asymmetric and symmetric at the side walls of the mixing channel. Each case has named as per the arrangement of hydro-active patches (hydrophobic/hydrophilic) in a micro-mixer. They were tested for various values of species decay constants allowed in the mathematical model. The flow is laminar at a Reynolds number of 0.5. Two meshes were employed both meeting the criteria of grid convergence tests for 5% error. It was found using comparison with analytical results that a finer grid shows spurious numerical features in the symmetric case. Surprisingly, this error was absent in the coarser grid. GCI test of simulations would thus not be sufficient for the purpose of validating the engineering design of micro-mixers. Sensitivity of the results to difference equation solving techniques has been addressed by setting suitable data structures (solvers and preconditioners) in OpenFOAM.
Apart from the main focus of our study, the developed model will be helpful to study mixing of the species having half-life and as a designing aid for the devices using decaying phases in fluid flow. It can be used to study the concentration of reactants at different time scales in nanoparticle synthesis [37]. It will also be helpful in particular engineering applications like ones involving microcapsules. Release of the active principles contained in microcapsules and its behavior over time into the medium such as in [38] can be simulated using our OpenFOAM model. The paper is organized as follows. Section II provides details of device design followed by the mathematical model described in section III. Numerical simulation is detailed in section IV leading to results on error and their discussion in section V. We present analytical proofs for error analysis in section VI and conclude in section VII.

II. DEVICE DESIGN
A passive micro-mixer with a rectangular cross section with hydro-active patches on the walls has been considered as a typical device design, for which simulations are carried out © 2019, IJSRPAS All Rights Reserved 12 and studied for numerical error. Figure 1 shows the schematic illustration. These three cases of patch arrangements are considered for the study which are termed as symmetric, asymmetric and no slip. They are named as per the arrangement of alternate hydrophobic and hydrophilic patches in a micro-mixer. The base geometry is rectangular with an aspect ratio of 4:1((cross-flow width is 1000 micron and in-flow length is 4000micron). Two inlets are provided for two different fluids. The issues of the numerical errors are studied by changing the method of solving the difference equations involved. Boundary conditions related to specific design factors as well as mesh refinements have also been changed to address the impact of numerical error on reliability of simulation results.

III. GOVERNING EQUATIONS AND MATHEMATICAL MODEL
A mathematical model to analyze the fluid flow and the concentration profile, in the micro channel is presented here. The model is based on continuum flow. For the continuum flow regime, the familiar Navier-Stokes-Fourier (NSF) continuum-fluid equations provide an excellent approximation for fluid flows that are very close to equilibrium. It is assumed that local macro-properties can be described as averages over elements that are large compared to the microscopic structure of the fluid assuming a current density for amount of active ingredient as c D v c   (2) and invoking conservation of the active ingredient.
c is the concentration of one of the phases, v is the velocity field, D is the diffusion constant and ƛ is decay constant. Further, we assume incompressible flow for varying both velocity and concentration .The momentum balance for such a flow can be expressed as Where p is the pressure field, ρ is the density and η is the dynamic coefficient of viscosity. We have assumed that Navier Stokes equations hold in case of our two phase fluid.

IV. NUMERICAL SIMULATION
Simulation solver is developed using OpenFOAM C++ library directories viz. case, solver and utility by suitably engineering data structures contained in them. The information regarding the discretization (mesh), the values of parameters of partial differential equations (PDEs) and the side conditions from the model were implemented in the case directory. The flow fields and the PDE along with the nature of side conditions are declared using data structures in the solver directory. These were implemented using the mathematical model and customized for the application. The numerical method used for discretizing the partial differential equation was implemented in the case directory. Similarly data structures were set up concerning the numerical technique for solving the resulting difference equations which addresses the main issue in our investigation. The numerical methods the difference equation techniques are described in the subsequent sections. Gauss linear and Gauss upwind schemes are employed for gradient and divergence discretization in PDEs of the mathematical model. These appear to be sufficiently reliable to provide broad trends in the flow fields. Such as upwind schemes avoid spurious oscillatory solutions and Gauss linear schemes have errors which reduce as squares of cell separation [29]. We are interested in demonstrating Laplacian operators in the PDEs which are truncated using the Gauss linear corrected scheme for the similar reasons and likewise for the time derivative scheme which was chosen to by Euler.
Further the resulting difference equations were numerically addressed using pre-conditioned bi-conjugate gradient and bi-conjugate gradient stabilized techniques. The pre-conditioners such as Diagonal based Incomplete Cholesky (DIC), diagonal incomplete LU (DILU), and diagonal were used for pressure, concentration and velocity. The OpenFOAM solvers are iterative, every iteration produces a residual for each equation, and lower residual values implying better accuracy. This ensures proper iterative convergence given the matrix form of the equations using smoothers available for pressure and concentration.
Smooth solver technique was employed for the velocity difference equation and a Symmetric Gauss Seidel smoother was used. Our objective is to study the effect of stabilization of numerical technique as well as the effect of grid size on the predictions. We also seek to relate this to boundary conditions used in the PDE model. The rectangular geometry was used as its structured mesh was easy to create. Also use of rectangular structured mesh minimized the computational time. Structured mesh has high accuracy, the points could be accessed in all directions, and refinement was much faster, also suitable for multi-mesh. The post processing of the results has been relatively easy task because the logical arrangement of the grid coordinates made it straightforward to plot the results [39]. Mixing was completely dependent on the diffusion of the species at the interface between the two liquids. Such interface was available in rectangular type of geometry [31]. This mixing portion of the channel was divided lengthwise into four equal parts. The top and bottom surface of each part was divided into patch of equal length. Each patch can carry either a no slip or slip boundary condition for velocity. The side walls bore no slip conditions. Inlet velocity conditions were kept constant and equal volume flow rate at both inlets. Zero pressure at the outlet. Three cases with different combinations of the above conditions are simulated to study the effect on mixing performance.
We realized that the round off error and tolerance affect the performance of the solver. This is not captured in Grid Convergence Index based methods. We changed the difference equation solving technique to a more stable one and found that the error could be removed. We noted that the boundary conditions also affect such type of error which makes this issue crucial from the design point of view [28]. The structured meshes used for the simulation are given in figure 2.   Simulations of all the cases were run using the developed solver with decay constant for unsteady flow. The concentration equation was solved using PBiCG (preconditioned bi-conjugate gradient) solver and preconditioned using DILU (diagonal incomplete LU) for various decay constants. Concentration values obtained during the simulation were exported and plotted to study the mixing performance of the micro-mixer. Plots presented in figure 3 (a) shows the glitches for the decay constant (ƛ) of 0.14 and 0.142 while flat portion of the graphs indicate the uniform mixing in a micro-mixer. Perturbations due to hydro-active patches has been energized the flow which resulted enhancement in the mixing. We have shown in section VI that the glitches obtained are due to numerical error. We therefore solved the same equations using PBiCGstab (preconditioned bi-conjugate gradient stabilized) instead of using PBiCG and preconditioned using DILU (diagonal incomplete LU) for various decay constants shown in Figure 3 (b). No glitches are evident in the plots and the numerical error has been eliminated.
During the study it was found that PBiCGStab solver is faster and more stable than PBiCG as it uses two residual vectors. It is also more robust with DILU preconditioner. The preconditioner is also changed from DILU to diagonal and obtained results have been presented in Figure 2 (c). We find that the results match with Figure 2 (b) which confirms that choice of iterative solver is more important than that of preconditioner. To check the results and estimated errors obtained in symmetric case, we nonetheless have obtained the Grid Convergence Index as prescribed by Roache [32] , [33] using two grid levels and concentration as the dependent variable. We have one grid coarser than the other by a factor of 1.33, and we obtained an apparent order p=4.3, leading to an estimate of 5% accuracy in the concentration values. This compares well with the observed 3% variation of concentration variance and further supports the correction in the error. The grid convergence index was proposed as a diagnostic on numerical accuracy by Roache in [32].
(a) solver: PBiCG, preconditioner: DILU, (b) Solver: PBiCGStab, Preconditioner: DILU, (c) Solver: PBiCGStab, Preconditioner: diagonal , No anomalies like glitches are observed by coarsening the grid for the same boundary conditions. Graphs in each frame are identical up to scaling as theoretically expected  Figure 4 shows the results of symmetric case with the mesh cell size of (100*100*100) shown in Figure 2 (b). It was found that the glitches were eliminated even for the decay constants of 0.14 and 0.142 in the case a where the concentration equations were solved using PBiCGsolver (preconditioned bi-conjugate gradient) and preconditioned using DILU (diagonal incomplete LU). The problem occurred for the finer grid which suggests that considerable numerical © 2019, IJSRPAS All Rights Reserved 15 error may be caused due to cell sizes and cell alignment in the mesh. The nature of the graphs in Figure 4 is not smooth as compared with the finer grid (graphs in Figure 3) but we found the same mixing trend. Occurrence of the error in the finer grid shows that only carrying out tests using Grid Convergence Index would not be sufficient for checking the suitability of such simulations as design aids.

B. Case I: Asymmetric
(a) solver: PBiCG, preconditioner: DILU, (b) Solver: PBicGStab, Preconditioner: DILU, (c) Solver: PBicGStab, Preconditioner: diagonal, No glitches are observed for these boundary conditions.   Figure 1 (b). The simulations of all cases were run using the same numerical methods as in the symmetric case. The glitches were not seen in this case but mixing is not prominent as compared with symmetric boundary conditions. Still as compared to no slip boundary conditions, mixing is more as the flat portion can be realized in the center of the graph. It proves that perturbations due to the alternate hydroactive patches are stimulating the two species to mix together. We did not find any numerical error in the graphs in this case therefore the simulations for the other grid were not carried out. This again suggests that boundary conditions too are affecting the numerical error.  This example highlights the absence of hydrophilic and hydrophobic patches in the micro-mixer as in Figure 1 (c). It is a plain micro-mixer which is used to compare the results of the mixing due to the addition of the hydro-active patches (geometry) which causes the agitation in the flow. Again all the simulations were run using the same combination of numerical method and preconditioners. There is no flat portion in these graphs which indicates very less mixing during the flow. In this case also we did not see any glitches in the graphs so we did not carry out the simulations for the other grid. Absence of the glitches shows the importance of considering together the effects of numerical method, grid size and boundary conditions during simulations. The numerical error observed to occur only in the symmetric boundary condition case which has been seen to offer better mixing performance if used in design of micro-mixers. Further, the error disappears upon coarsening the grid in contrast to refining the grid as suggested by Roache criterion (at anticipated 5% error). We believe that it could be because of round off error as the unstabilized difference equation algorithm requires more iterations, a situation worsened if more cells exist in the mesh. Indeed, as we found, stabilization eliminates the error for the finer mesh.

VI. ERROR ESTIMATES: REQUIRED ANALYTICAL RESULTS
We apply an integrating factor of t e  to the Equation (1) to which is the PDE model for concentration c when the decay constant vanishes. A scale factor depending on decay constant and time can be multiplied to obtain the solution of all concentration profiles. This shows that all the graphs should have the same shape under scaling. No glitches are expected as seen in the symmetric case for the fine grid and use of PBCiG solver if there had been no numerical error.

VII. CONCLUSION AND FUTURE SCOPE
Numerical errors in simulations of passive micro-mixers are important to consider for predicted mixing performances of device designs. We consider a design of mixing channels possessing hydro-active patches on the top and bottom floor. Simulations for asymmetric, symmetric arrangement of the patches and for plain wall conditions were carried out using OpenFOAM CFD software. The simulations show features (glitches) in spatial concentration profiles for certain conditions. Based upon analytical treatment we establish that these features are due to numerical error and are therefore spurious. These spurious features are observed even though the mesh employed for the above simulations meets the criteria of grid convergence test for 5% error. Surprisingly, the spurious features were absent in simulations carried out using a coarser grid. GCI tests of simulations would thus not be sufficient to assure reliability of such simulations for the larger purpose of validating the engineering design of micromixers.
Further, we find that errors disappear not only after changing mesh sizes but also upon changing the numerical technique employed for solving difference equations suggesting a strong interplay between these factors. A possibly important role played by round off errors is thus indicated by our study. This effect is also sensitive to the kind of boundary conditions imposed, which makes the above interplay very crucial for design purposes. Instead of Grid Convergence Index, we suggest that a better criterion is required from the design point of view. Such a criterion should be additionally informed by difference equation solving techniques and boundary conditions in microflow simulations.
These passive micro-mixing simulation cases could be used for benchmarking, but we note here that they would also be useful to study particular micromixing applications.
Chemical decay of species in unsteady micro-flows is directly addressed by this model and its effects can be studied when a mixing ingredient has half-life comparable to flow time scales in the device.