Uncertainty Analysis of Mixing Efficiency Variation in Passive Micromixers due to Geometric Tolerances

The geometric layout is the key factor for enhancing the efficiency of the fluidmixing in passivemicromixers.Therefore, by adjusting the geometric design and by controlling the geometric parameters, one can enhance the mixing process. However, through any fabrication process, the geometric parameters present slight, inherent variation from the designed values than might affect the performance of the micromixer. This paper proposes a numerical study on the influence of the unavoidable geometric tolerances on the mixing efficiency in passive micromixers. A probabilistic simulationmodel, based on theMonte Carlo method, is developed and implemented for this purpose. An uncertainty simulationmodel shows that significant deviations from the deterministic design can appear due to small variations in the geometric parameters values and demonstrates how a more realistic mixing performance can be estimated.


Introduction
One of the challenges in microfluidics is the mixing of fluids.Due to the low Reynolds number flows, the mixing process is done mainly by diffusion, which makes it very slow.In order to achieve a uniform concentration by mixing two fluids in a microchannel with the width of the order of hundreds of micrometers, a length of order of centimetres is needed [1].Numerous methods for enhancing the mixing process by shortening the mixing length and the mixing time have been proposed in the literature, among which can be mentioned the geometry effect, hydrodynamic focusing, alternate injection, electrokinetic method, droplet mixing, and stirring by particles [2].The most simple approach is the one based on the geometric effect.The mixing length can be shortened by modifying the geometry of the microfluidic channel.Several geometric designs for mixing enhancement have been proposed in the literature: microchannels with serpentine elements [2,3], two-layered crossing channels [2,4], and microchannels with grooves or obstacles on the bottom wall [2,[5][6][7].
By increasing the geometric complexity, also the complexity and difficulty of the fabrication process are being increased.A more complex fabrication process implies more variations of geometric parameters due to the fabrication tolerances.Most of the techniques developed up to date for monitoring the effects of the variation of geometrical parameters due to the fabrication tolerances on the behaviour of microfluidic systems are based on experimental methods [8,9].The experimental testing of the device is not always convenient, since many applications like biomedicine, pharmaceutics, or environmental monitoring require the development of sterile or disposable devices and physical testing would imply the device contamination.On the other hand, for some applications there are no experimental testing procedures available.For example, for mixing microfluidic devices there are no standard methods for determining the uniformity of concentration across a microchannel, in order to assess the quality of the mixing process.
In other applications like biophysics or biochemistry, functionalised surfaces are needed and the experimental testing would imply the damaging of these surfaces.Another approach has been to integrate testing mechanisms in the device (like optical sensors, capacitance measurement electrodes, impedance sensors, etc.).This approach brings along two main disadvantages: the design complexity and implicitly the cost of the device are increased and, on the other hand, there are needed supplementary fabrication steps, which are not always compatible with the materials used for manufacturing the microfluidic device.
Given the unique challenges presented by the experimental testing of microfluidic devices, in this work, the development of a simulation based methodology for estimating the influence of geometrical parametric deviations is proposed, which is caused by the fabrication tolerances, on the performance of passive micromixers.Passive micromixers are devices that use the geometry as the mixing enhancement factor.These types of devices are cost effective, from the fabrication point of view, and also easy to integrate in complex microfluidic systems; therefore, they have a high potential of applicability.The purpose of the study presented in this paper is to show how inherent geometric uncertainties generated by the fabrication processes propagate through the design and influence the final performance of the device.A demonstration of the proposed simulation methodology is done using as test structure a passive micromixer with modified tesla elements [10].

Geometry and Description of the Passive Micromixer Test Structure
The layout of a micromixer with a modified tesla element is presented in Figure 1 [10].The fluids are introduced in the device through two separate inlet channels that meet in a main microchannel.A secondary channel leaves the main channel at an acute angle, while the main channel has a 180 ∘ turn that forms the tesla element and meets again with the secondary channel.The fluids that are introduced in the micromixer are water, through the upper inlet, and a NaCl solution with a concentration of 9 mg/mL (154 mol/m 3 ), through the lower inlet.The geometry of a tesla element is characterised by the following geometric parameters (Figure 2): width of the main channel ( 1 = 500 m), height of the main channel (ℎ 1 = 500 m), radius of the inner curve ( = 125 m), opening of the secondary channel (sc = 100 m), width of the outlet channel ( 2 = 500 m), and height of the outlet channel (ℎ 2 = 500 m).The widths of the two inlet channels are given by the parameters  1 =  2 = 500 m.

Numerical Simulations
The deterministic modelling of the mixing process inside the micromixer with tesla elements embeds a fluid dynamics model for the flow of the fluids through the micromixer and a diffusion process for the spread of molecules through random motion from regions of higher concentration to regions of lower concentration.The fluid flow through the micromixer is modelled with the Navier-Stokes model for incompressible fluids which describes the flow of viscous fluids with momentum balances for each component of the momentum vector in all spatial dimensions, under the assumption that the density and viscosity of the modelled fluid are constant.
The Navier-Stokes model that describes the incompressible fluid flow model is given by the following system of equations [11]: where  is the fluid's density (kg/m 3 ),  represents the velocity vector (m/s),  is the pressure (Pa),  is the dynamic viscosity (Pa⋅s), and  is the identity matrix.The boundary conditions that have been considered for the Navier Stokes model are zero velocity at the wall interface, an initial velocity of 10 m/s at the two inlets, and zero pressure at the outlet.The diffusion process is modelled by Fick's law that describes the diffusive transport in the flux vector.The convection and diffusion equation is [12] where   is the concentration of species  (mol/m 3 ),   denotes its diffusion coefficient (m 2 /s), and  refers to the velocity (m/s).The velocity is obtained from the Navier-Stokes model for incompressible fluids.The boundary conditions for the diffusion model are the given concentrations of the fluids at the inlets: 0 mol/m 3 , respectively, 154 mol/m 3 and the convective flux condition at the outlet, which means that mass is transported out of the domain by convection only.The deterministic simulation of the passive micromixer with modified tesla elements has been performed with the COMSOL simulation software [13,14].In order to apply both models needed for the description of the mixing process a coupled simulation is needed.The coupled simulation is performed by using the Incompressible Navier-Stokes application mode of the MEMS module of COMSOL and the Convection and Diffusion application mode of the Multiphysics module.Because the geometry of the micromixer is prismatic (it has a constant section along the 0 axis), a 2D simulation is sufficient for the study of the mixing process.A summary of simulation settings and input parameters used for performing the numerical simulation of the mixing process inside the micromixer with tesla elements is presented in Table 1.
First a simulation for a micromixer structure with one tesla element has been performed.The results of the Incompressible Navier-Stokes application mode show an increased velocity at the intersection of the U-turn of each tesla element with the secondary channel (Figure 3).The streamline plot and the arrow graph of the velocity field show that the fluid tends to follow the inner curved surface.It can be seen that in the secondary channel there are two streams of fluid that meet from opposite directions.The resulting values of the -velocity and the -velocity along the surface of the micromixer are used as input data for the Convection and Diffusion application mode.This second application mode calculates the concentration distribution along the 2D profile of the micromixer.A plot of the concentration distribution along the micromixer overlapped by a streamline and arrow graph of the velocity field (Figure 4) shows that the stream of fluid with 0 mol/m 3 concentration comes out of the U-turn of the tesla element with a higher velocity than the stream of fluid with the 154 mol/m 3 concentration that is coming from the secondary channel.These two streams of fluid come from opposite directions and meet in the secondary channel.In the area where the two streams of fluid meet, the diffusion process is accelerated.
The concentration distribution graphic from Figure 4 shows that one tesla element is insufficient for efficient  mixing.In order to obtain a more uniform concentration in the outlet channel, it is necessary to increase the number of tesla elements.The tesla elements are connected in series such that the outlet of one tesla element coincides with the inlet of the following tesla element.In order to determine the minimum number of tesla elements required for an efficient mixing, further analyses have been performed.Further results are presented for micromixers with 8 connected modified tesla elements.The concentration distribution resulting from the Convection and Diffusion application mode of COMSOL for a micromixer with 8 connected tesla elements is presented in Figure 5 [13].It can be observed that along the seventh and the eighth tesla elements of the micromixer the concentration is visibly more uniform than in the previous tesla elements.

Mixing Efficiency
In order to analyse the evolution of the concentration distribution along the micromixer with 8 tesla elements, a graph with the cross section plots determined in the outlet channel of each tesla element is realised (Figure 6).It can be observed that the concentration distribution varies significantly between the two walls of the micromixer at the outlet of the first four tesla elements, while the cross section profiles of the concentration distribution at the outlet channel of the sixth, seventh, and eighth elements do not present significant variations.From this it could be concluded that a design of the passive micromixer with six tesla elements should ensure a sufficient uniformity of the concentration distribution across the channel.However, in order to have a more rigorous  decision on the most efficient design, a mixing efficiency coefficient needs to be defined.
The mixing efficiency is evaluated according to the following formula [15,16]: where  represents the mixing efficiency,  is the concentration distribution across a transversal section of the micromixer,   represents the concentration of complete mixing,   represents the initial concentration, calculated at a cross section of the micromixer, near the inlets, and  represents the width of the considered cross section.The numerical values of the concentration distribution obtained by the finite element simulation are exported in Matlab [17], and the mixing efficiency is computed using numerical integration.The results of the computed mixing efficiencies for all the tesla elements of the passive micromixer are presented in Table 2.It can be observed that, for the first tesla element, the mixing efficiency is very low, 36.65%, and that it is increasing with the number of tesla elements.Also, it can be observed that there is a very slight improvement of the mixing efficiency from the fifth up to the eighth tesla elements: from 95.98% it rises to 98.83%.These results show that a design with minimum five tesla elements should ensure a satisfactory mixing efficiency (above 95%).

Uncertainty Analysis of Mixing Efficiency Variation
In passive micromixers such as the tesla micromixer presented in the previous section, the geometry of the device is the main factor that contributes to the acceleration of the mixing process.However, in practical applications, the designed geometry suffers small variations due to the fabrication tolerances.The dimensions of the various geometric features used in the design can be obtained in manufactured structures only up to a certain precision.Small variations are inherent.
To verify how the fabrication tolerance influences the functionality of the passive micromixer with tesla elements, a Monte Carlo simulation is performed.The purpose is to study how the variations in the geometry propagate through the model and affect the final result.The parameter analysed in the Monte Carlo analysis is the mixing efficiency coefficient (3).
The principle behind the proposed probabilistic simulation method is that the geometrical parameters are no longer considered to be constants, but they are implemented as random variables that follow given probability distributions.For the implementation of the geometry variation caused by the fabrication tolerances, it is necessary to parameterize the geometry of the studied microfluidic system.Then, the geometric parameters are modelled using probability distributions.Once the variation laws of the geometrical parameters are defined, the Monte Carlo simulation can be applied [18,19].For one iteration of the Monte Carlo simulation, a random set of input geometrical parameters is generated from their probability distributions and the numerical simulation is carried on using these input values.The result of the numerical simulation is stored in a collection of output data.These steps are repeated until a given number of iterations has been performed.In the end, the output data is analysed by statistical methods.

Geometric Variations.
Before applying the Monte Carlo simulation, parameterization of the geometry of the micromixer is done.The purpose of this study is not to focus on a particular fabrication process, since some geometric tolerances are present regardless of the chosen technology and materials.Therefore, it has been considered a normal distributed variation of up to 5%, for all the geometric parameters involved in the design of the considered micromixer.
Figure 2 presents all the geometric parameters necessary for the description of one tesla element.Each tesla element is characterised by three geometric parameters: the width  of the main channel (), the radius of the inner curve () and the opening width of the secondary channel (sc).For the parameterization of the micromixer with 8 connected tesla elements, there are necessary 24 parameters for the tesla elements and in addition two parameters for the widths of the inlet channels and one parameter for the width of the outlet channel, resulting in a total of 27 parameters.Table 3 presents the list of the 27 parameters used for the definition of the micromixer geometry with 8 tesla elements and the law of variation for each parameter.All the parameters variations follow a normal distribution, with the standard deviation, , calculated such that the variation range is ±5%.

Monte Carlo Simulation.
The Monte Carlo simulation is executed using this parameterization model of the geometry of the micromixer.For an iteration of the Monte Carlo algorithm, a set of random parameters is generated using the probability laws presented in Table 3.The simulation of the micromixer is executed like a deterministic finite element simulation having as entry data the geometry defined by the randomly generated set of parameters and the mixing efficiency coefficient is calculated.Several iterations are performed and afterwards the obtained data set is analysed statistically.
In order to have a more precise result of the statistical analysis it is preferable to have a volume of data as large as possible.Therefore, a trial of 100 iterations is performed.For testing the convergence of the Monte Carlo simulation an intermediate mean values plot is performed with the corresponding confidence interval (Figure 7).The convergence plot shows that no significant variations of the mean value of the data set are obtained after about 50 iterations.This means that the considered trial of 100 iterations of the Monte Carlo simulation is sufficient to provide an accurate statistical analysis of the results.

Results and Discussions.
The Monte Carlo simulation method has been used to compute the mixing efficiency parameters at the outlet channels of the final three tesla elements of the passive micromixer, for normal distributed random samples of the geometric parameters.Three sample data sets have been obtained for the four mixing efficiency coefficients.Figure 8 present the histogram representation of the distributions of these data sample sets.It can be observed that the distribution of the mixing efficiency parameters is not normal; therefore, for the statistical analyses of these data, nonparametric methods are required.
Since the objective of this study is to quantify the performance of the passive micromixer, regardless of the fabrication tolerances, further, there are computed prediction intervals for the mixing efficiency coefficients from the data samples obtained by Monte Carlo simulation.A nonparametric method, based on order statistics [20], is used for this purpose.The th order statistic of a statistical sample is equal to its th-smallest value.If  () and  () are order statistics of the sample with  <  and  +  =  + 1, then [ () ,  () ] is a prediction interval for  +1 with significance level equal to ( + 1 − 2)/( + 1).Hence, if a prediction interval with a significance level higher then  is sought, then the lower limit of the prediction interval can be found by choosing the th order statistic with  = int((1/2)(+1)(1−)).The results for computing the prediction interval for the three mixing efficiency coefficients are shown in Table 4.
The first line of Table 4 displays the value of the mixing efficiency found by deterministic simulation.The second line displays the lower boundary of the prediction interval and the third line displays the upper boundary.It can be observed that the value of the mixing efficiency can vary with up to two percent from the value computed by deterministic design.The prediction intervals offer a more realistic estimation of the expected mixing efficiency in a fabricated micromixer.
Using the data sampled obtained from the Monte Carlo simulation, the empiric cumulative distribution functions have been determined using Matlab [17] and they have been used to compute the probability that the mixing efficiency will fall under 95% in all the three considered cases.The results show that variations of up to 5% of the geometric parameters of the passive micromixer lead to a probability of 11.1% that the mixing efficiency will be lower than 95% for devices with a design with 6 tesla elements, while for a design with 7 tesla elements this probability is 2.08%.All these considered, it can be concluded that a design with at least 7 tesla elements would be more reliable than a design with 6 tesla elements.

Conclusions
In this paper, a study of the influence of geometric parameters uncertainties on the performance of the fluid mixing process in passive micromixers has been presented.Such a study is useful, because in these types of structures, the mixing process is influenced mainly by the geometric effect.Therefore, it is important to predict how inherent variations of the geometry caused by the fabrication tolerances may influence the behaviour of the device.The structure chosen as a demonstrator is a passive micromixer with modified tesla elements due to the complexity of its geometry and the high number of parameters needed for defining the geometry.The study shows that the Monte Carlo simulation offers supplementary information on the behaviour of the device to the deterministic finite element modelling.By statistical analyses of the data samples provided by the Monte Carlo simulation, it was possible to determine a range of variation of the mixing efficiency.The results also show that variations of up to 5% of the geometric parameters of a passive micromixer with 6 tesla elements lead to a probability of 11.1% that the mixing efficiency will be lower than 95%.One possible solution to overcome this issue would be to take into consideration a design with 7 tesla elements, for which the probability that the mixing efficiency would become lower than 95% is only 2.08%.

Figure 1 :
Figure 1: Layout of a passive micromixer with tesla elements.

Figure 4 :
Figure 4: Concentration plot overlapped with streamline plot and arrow graph of the velocity field for the micromixer with one modified tesla element (Comsol simulation).

Figure 5 :
Figure 5: Comsol simulation of the concentration distribution along a micromixer with 8 connected tesla elements [13].

Figure 6 :
Figure 6: Cross section concentration plots at the end of each of the 8 connected tesla elements of the micromixer.
the confidence interval Intermediate mean value Upper boundary of the confidence interval

Figure 7 :
Figure 7: Convergence plot of the Monte Carlo simulation for the micromixer with 8 connected tesla elements.

Figure 8 :
Figure 8: Histograms of the results of the Monte Carlo simulation of the tesla micromixer at the outlets of the 6th up to the 8th tesla elements.

Table 1 :
Comsol simulation settings and input parameters used for performing the numerical simulation of the mixing process inside the micromixer with tesla elements.

Table 2 :
Mixing efficiency coefficients computed at the cross sections of the outlet channels of every tesla element of the passive micromixer.

Table 3 :
Variation models of the geometrical parameters of the micromixer with 8 connected tesla elements (for each parameter,  is calculated such that the variation of the parameter is up to 10% from its nominal value).

Table 4 :
Prediction intervals for the mixing efficiency coefficients computed along transversal cross sections at the outlet channels of the final three tesla elements of the passive micromixer.