Computational Evaluation of Mixing Performance in 3-D Swirl-Generating Passive Micromixers

Abstract

Computational Fluid Dynamics (CFD) tools are used to investigate fluid flow and scalar mixing in micromixers where low molecular diffusivities yield advection dominant transport. In these applications, achieving a numerical solution is challenging. Numerical procedures used to overcome these difficulties may cause misevaluation of the mixing process. Evaluation of the mixing performance of these devices without appropriate analysis of the contribution of numerical diffusion yields over estimation of mixing performance. In this study, two- and four-inlet swirl-generating micromixers are examined for different mesh density, flow and molecular diffusivity scenarios. It is shown that mesh densities need to be high enough to reveal numerical diffusion errors in scalar transport simulations. Two-inlet micromixer design was found to produce higher numerical diffusion. In both micromixer configurations, when cell Peclet numbers were around 50 and 100 for Reynolds numbers 240 and 120, the numerical diffusion effects were tolerable. However, when large cell Peclet number scenarios were tested, it was found that the molecular diffusivity of the fluid is completely masked by false diffusion errors.

1. Introduction

The use of micromixers in microfluidic devices has become popular due to multiple advantages of small scales, less energy requirement, low cost, and high efficiency [1]. Micromixers are essential components of microfluidic devices [2] which are used in various fields for analysis, diagnostics, control, and monitoring [3]. The basic function of micromixers is to provide rapid mixing of two or more fluids over a short distance. These devices are characterized as active and passive micromixers [1,2,4,5] based on their external energy use. While active micromixers require an external disturbance (e.g., electrical, magnetic, sound etc.) and present high mixing efficiencies in short distances [3], passive micromixers function using the fluid flow energy within the micromixer and utilize two mixing phenomena—advection and diffusion—in relatively less complicated mixing structures [6]. In passive micromixers, achieving rapid mixing in short distances requires special geometric designs due to highly laminar flow characteristics (e.g., Re < 100 [6]) and small molecular diffusion coefficients on the order of 10⁻⁹–10⁻¹¹ m²/s [7]. At extremely low Reynolds numbers (e.g., Re < 1), mixing is mainly controlled by molecular diffusion [8] and obtaining good mixing requires long mixing channels and time. On the other hand, increasing the flow rate (e.g., Re > 1) creates advection dominant flow within the micromixer whereby mixing becomes even more challenging unless the contact surface area between fluid bodies is increased by creating secondary flows in the micromixer. In this case, high flow velocities along with low molecular diffusivities create complexity in terms of numerical evaluation of mixing using a Computational Fluid Dynamics (CFD) approach which is extensively used in passive micromixer studies [9,10,11].

In the numerical solution of advection dominant systems, numerical problems usually occur as oscillations on the concentration front and back (i.e., lack of dispersion for sharp concentration gradient flows) and unphysical diffusion (i.e., false diffusion or numerical diffusion). Most common numerical methods, such as Finite Volume Method (FVM) and Finite Element Method (FEM), that are employed in CFD applications suffer similar numerical error outcomes. As discussed in detail in Reference [11], these numerical techniques use different approaches to solve 3-D fluid flow and scalar transport equations to overcome numerical complications. However, these applications also contribute to false diffusion in the solution when Pe ≫ 1, which needs to be evaluated and reported.

In the current literature, only a few papers [9,10,11,12] have focused on the influence and extent of numerical problems in the estimation of mixing performance in passive micromixers. Other than these, artificial mixing effects of the numerical methods have been overlooked and erroneous mixing performance results have been reported in several studies. For example, Cortes–Quiroz et al. [13] numerically investigated fluid mixing in a swirl-generating 3-D passive T-mixer with a 2000 µm long mixing channel for a Re number range between 100 and 500 using the FVM. To examine the grid sensitivity at Re = 250, they employed six different mesh levels which used structured brick elements with mesh densities between 3 million and slightly over 5 million and selected a mesh around 4.3 million elements for all simulations in the study. Although very high Peclet numbers (e.g., on the order of 10⁵) were studied, numerical diffusion effects on mixing were not reported. Modifying the inlets and constricting the mixing channel, Zhang et al. [14] studied similar 3-D swirl-induced T-mixer for various Reynolds numbers between 10 and 70 using FEM [13]. All simulations were conducted using a mesh density around 277,000 with structured and unstructured elements. The numerical results of the original 3-D micro T-Mixer (OTM) in this study were compared with an experimental study [15]. While using a 5000 µm shorter mixing channel in numerical simulations, approximately 20% higher mixing efficiency was reported when compared with the experimental result at Re = 70. Such a large discrepancy was more likely to have emerged due to high numerical diffusion, which was due to the use of a rather coarse mesh in simulations.

Ansari et al. [16] surveyed a comparative analysis of classical T-mixer and vortex T-mixer both numerically and experimentally for Re numbers between 1 and 80. In this study, numerical simulations were conducted using a structured hexahedral mesh with a grid size of 4 µm at the T junction and larger elements were used at the downstream end with a total element number of approximately 1.3 million. In this numerical solution, FVM was used. Even though it was realized that numerical simulations were not free of numerical diffusion, this was not examined and reported further.

Izadpanah et al. [17] studied fluid mixing for Re cases 75–400 in a T-shaped and double T-shaped micromixer geometry and investigated the effects of vortex and engulfment flow regimes on mixing using FVM. For a mesh independence study at Re = 260, four different mesh densities consisting of hexahedral elements ranging from 0.24 million to 1.3 million were tested and a mesh density close to 1 million elements was selected for simulations. At such a low mesh density and high Re cases, numerical solutions will be exposed to high false diffusion, but authors evaluated the mixing performance without mentioning the contribution of numerical diffusion effects.

Chen et al. [18] investigated micro-scale mixing using swirl-generating designs for several inlet and mixing box configurations up to Re = 100 employing FVM and structured hexahedral mesh. The authors mentioned that to reduce numerical diffusion effects, mesh refinement was maintained until the mixing index difference was below 5% between mesh levels, but they did not report the employed grid sizes and mesh density levels in the simulations.

In Reference [11], the numerical diffusion effects of FVM were exhaustively investigated in a simple T-shaped passive micromixer and it was found that maintaining a good mesh flow alignment is highly important to reduce false diffusion in advection dominant transport systems. This was only possible by using structured hexahedral elements in the solution of the transport equation. In contrast, when prism and tetrahedral elements were used, physical diffusion in the micromixer was severely masked by numerical diffusion errors which manifested themselves as unphysical mixing in numerical simulations. Unlike the simple stratified and vortex flow regimes observed in simple T-shaped passive micromixers in References [11,19,20], swirl flows are prone to create high numerical diffusion effects as a result of the continuously changing flow direction of fluid pairs. Even if a structured hexahedral grid type is used for swirl flow systems, constant disorientation of flow and cell boundaries will lead to high amounts of false diffusion, which contributes to numerical mixing. Hence, the numerical solutions of these systems require extreme caution to avoid numerical errors which may alter the interpretation of the physical mixing problem entirely.

Thus, it is essential to characterize the swirling passive micromixers in terms of the limits of numerical diffusion errors in CFD simulations to provide reliable mixing evaluations. In this study, two different 3-D swirl-generating passive micromixer designs with two- and four-inlet constructions were investigated for various fluid flow and transport scenarios using the FVM, and the false diffusion generated in numerical solutions was quantified to show the adverse effects on mixing performance estimates. It should be noted that the extent of the present research is not only limited to swirling flow systems, but also for advection dominant micro-scale mixing systems in which there is a continuous flow direction change.

2. Governing Equations

The CFD tool is commonly used to mathematically describe fluid flow and transport in numerical passive micromixer studies. Laminar fluid flow and scalar transport equations for incompressible, isothermal, and conservative fluids under negligible gravitational effects can be defined as given in Equations (1)–(3) in which $\overrightarrow{U}$, ρ, p, µ, c, and D denotes velocity vector (m/s), density of fluid (kg/m³), pressure (Pa), dynamic viscosity of the fluid (Pa·s), solute concentration (mol/m³), and molecular diffusivity (m²/s), respectively.

$$\nabla ·\overrightarrow{U}=0$$

(1)

$$\rho \left[\frac{\partial \overrightarrow{U}}{\partial t}+\overrightarrow{U}·\nabla \overrightarrow{U}\right]=-\nabla p+\mu {\nabla }^2\overrightarrow{U}$$

(2)

$$\frac{\partial c}{\partial t}+\overrightarrow{U}·\nabla c=D{\nabla }^{2}c$$

(3)

In this study, mixing on a plane is calculated using the mixing index (MI) [14,21,22] as given in Equation (4). MI has a range between 1 and 0 which represents full mixed and unmixed conditions respectively. It should be noted that the MI approach has been employed due to uniform scalar distribution and symmetrical flow profiles observed on the planes across the streamwise direction. In non-uniform cases, the flow average mixing method discussed in Reference [19] should be used for accurate quantification of the mixing efficiency.

$$\mathrm{MI}=1-\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(\frac{c_i-c_{avr}}{c_{avr}}\right)}^2}}$$

(4)

where, c_i is the concentration at sample point i, c_avr is the average concentration on the given plane, and N is the total of sample points on the plane. Reynolds (Re) number and Peclet (Pe) number, given in Equations (5) and (6) respectively, are used to describe flow and transport characteristics in the micromixer respectively.

$$\mathrm{Re}=\frac{u{D}_h}{\nu }$$

(5)

$$\mathrm{Pe}=\frac{u{D}_h}{D}$$

(6)

where, u, ν, and D_h represent average flow velocity (m/s), kinematic viscosity (m²/s), and characteristic length (m) respectively. Characteristic length is defined as D_h = 2HW/(H + W) in which H and W are the channel height and width of a rectangular duct.

In the present study, the numerical solution of fluid flow and scalar transport in the micromixer is obtained employing OpenFOAM [23] (v5.0, OpenFOAM Foundation, OpenCFD Ltd., Bracknell, UK, 2015), which is an FVM-based solver. To solve the steady-state laminar flow and scalar transport equations, built-in solvers simpleFOAM and scalarTransportFOAM in OpenFOAM library were utilized respectively. In both solvers, advection and diffusion terms were discretized using second-order accurate upwind [24] and second-order accurate central difference schemes respectively. Simulations were terminated when the difference between the residuals fell below 1 × 10⁻⁶, which was assumed to be converged as usually practiced in several studies [16,17,19,25]. All simulations were conducted using the same computer configuration as in Reference [11].

3. Numerical Problem Description

The origin of numerical problems in passive micromixer studies is in the advection dominant transport applications (i.e., Pe ≫ 1) as extensively discussed in Reference [26]. Similarly, momentum and thermal energy transport may also present numerical difficulties. However, relatively higher kinematic viscosities and thermal diffusivities [9] reduce the numerical errors significantly for such applications.

In FVM, numerical diffusion develops in the discretization of the advection term of the transport equation [26]. The magnitude of this error changes depending on the accuracy of the numerical solution scheme employed, flow velocity, grid size, and the angle between flow velocity and grid boundaries [9,10,11]. In Figure 1, the extent of false diffusion effects in FVM are qualitatively shown for various discretization schemes, grid sizes, and mesh-flow alignment scenarios. In the literature, this is also known as the standard test as described and used in Reference [9]. In a 2-D square domain with an edge size of 1000 µm, a constant flow field was set across the domain for three different cases and each case was also tested for two different mesh levels, L1 (400 elements) and L2 (10,000 elements), using numerical schemes of first- and second-order accurate upwind cases. While Figure 1a,b show 45° and 90° flow-grid alignment for equally-spaced square grid type respectively, Figure 1c shows randomly changing flow orientation for a triangular mesh system. The constant flow fields applied for scenario (b) and scenarios (a) and (c) are (u, v) = (0.01 m/s, 0) and (0.01 m/s, 0.01 m/s) respectively, which corresponds to Re = 10 as calculated in Reference [9]. Thus, Figure 1 shows a steady-state scalar transport solution of a pure advection system (i.e., D = 0, Pe → ∞) using Equation (3) with imposed scalar values 0 and 1 used in inflow boundaries as shown in Figure 1a–c. The gradient of scalar value at all other boundaries were set to zero.

As can be seen from Figure 1b, numerical diffusion does not exist when fluid flow is orthogonal to the grid boundaries even if the first-order accurate upwind scheme and coarser grid elements are used. In contrast, when flow is oblique to the grid lines, the solution creates numerical diffusion depending on the angle between flow and grid lines. As shown in Figure 1a, when this angle is constant at 45°, the numerical solution creates the maximum amount of false diffusion on the coarser grid, L1. When triangular elements are employed, the produced false diffusion is less than the worst-case scenario (a) at L1, but the solution still contains a substantial amount of false diffusion as a result of randomly changing flow and grid alignment as shown in Figure 1c. Thus, less numerical diffusion with the use of triangular elements is due to relatively less obliqueness between flow direction and element boundary as a result of random orientation of the boundaries which may yield orthogonality for some elements throughout the domain. On the other hand, it is notable that using a second-order accurate scheme for discretization of the advection term and grid refinement significantly reduces the amount of numerical diffusion generated for scenarios (a) and (c). It should be pointed out that using a very fine mesh in passive micromixer studies may not be practical and possible due to increasing computational costs. The same 2-D standard test was extended to a four-inlet case as shown in Figure 2. In this case, the numerical diffusion generated is increased as a result of enlarged contact surfaces between fluid bodies. It can be seen in Figure 2d that the first-order accurate scheme along with a coarse mesh yielded a completely different result than the physical problem as a result of the constantly maintained 45° grid-flow angle inside the computational domain and increased fluid pairs in the transport domain. On the other hand, when Figure 2d is compared with Figure 2f,g, the use of triangular elements resulted in less false diffusion because of reduced inclination between the flow and grid boundaries. Moreover, applying a unidirectional flow slightly diminished the false diffusion production as shown in Figure 2g. This occurred because of improved mesh-flow orthogonality within the transport domain. The four-inlet numerical simulations produced higher amounts of numerical errors than the two-inlet solutions under the same flow conditions in this problem. The use of the second-order discretization scheme and smaller grid elements significantly reduced the numerical errors in scenarios (d), (f), and (g).

Therefore, numerical diffusion generated in numerical simulations of scalar transport depends on several parameters as shown the in 2-D standard test and all these factors need to be considered together to control and evaluate the amount of numerical diffusion generated. The 3-D systems create relatively complicated numerical diffusion patterns due to the contribution of additional dimensions and secondary flows. Although grid refinement may be a unique solution to eliminate numerical diffusion errors, since maintaining a constant mesh-flow orthogonality is not possible for most problem types, drastically increased computational requirement often prevents this approach. Thus, quantification of numerical diffusion becomes essential to analyze and report the mixing performance results properly.

4. Micromixer Design

To investigate the numerical diffusion effects in swirl-generating micromixers, two similar unobstructed 3-D passive micromixers were designed with two- and four-inlet configurations as shown in Figure 3. Micromixers consist of three branches as follows. A square mixing channel with an edge size of 100 µm and length (L) of 1900 µm, a 400 µm-wide square mixing box with a 100 µm height, and 500 µm-long inlet channels with a cross-section of 100 µm × 100 µm. The physical properties of mixing fluids are identical with a density (ρ) of 1 × 10³ kg/m³ and viscosity (µ) of 1 × 10⁻³ Pa·s. Boundary conditions used to obtain a steady flow field in simulations are zero pressure at the outlet, uniform velocity at the inlets, and no slip condition on the walls of the micromixers. In all scenarios, an equal amount of fluid is injected from the inlets which depends on the Reynolds number. Solute transport simulations are investigated by applying molar concentrations, 0 and 1 mol/m³, from the inlets as shown in Figure 3 with blue and red arrows respectively. For all cases, the first stationary flow field was obtained by solving Equations (1) and (2) simultaneously and later the flow field obtained was used in scalar transport simulations using Equation (3).

The computational domains in both the four-inlet and two-inlet designs were meshed using structured hexahedron-type elements at six different mesh density levels as given in

Table 1. Mesh properties and densities for six different levels.

Mesh Level	Δx (µm)	Total Mesh Elements (Million)
		Four-inlet	Two-inlet
M1	1.25–2.50	16.40	15.10
M2	1.50–2.50	10.00	9.20
M3	1.75–3.00	6.30	5.70
M4	2.00–3.00	4.60	4.10
M5	2.25–3.50	3.00	2.75
M6	2.50–3.50	2.50	2.10

. As shown in Figure 4, while smaller grid sizes in the second column of Table 1 were used in the mixing channel and partially in the mixing box, since high swirl generation is expected to occur in these sections of the micromixers, larger grid sizes were positioned in the inlet channels due to the fact that flows in the inlet channels are mainly unidirectional and mesh-flow alignment is kept seamlessly in these regions. To preserve a good mesh quality in the computational domain, the maximum aspect ratio within the hexahedron mesh elements was selected as 2.

To quantify and characterize the upper limits of the numerical diffusion effects, two different flow scenarios were tested with three molecular diffusion constants as given in

Table 2. Test cases for two-inlet and four-inlet micromixer designs.

Reynolds	240			120
DM (m2/s)	3 × 10−10	3 × 10−9	6 × 10−8	3 × 10−10	3 × 10−9	1.5 × 10−8
Peclet	8 × 105	8 × 104	4 × 103	4 × 105	4 × 104	2 × 103
Mesh Level	† PeΔ			† PeΔ
M1	10,000	1000	50	5000	500	100
M2	12,000	1200	60	6000	600	120
M3	14,000	1400	70	7000	700	140
M4	16,000	1600	80	8000	800	160
M5	18,000	1800	90	9000	900	180
M6	20,000	2000	100	10,000	1000	200

^† Average cell Peclet numbers are calculated using average velocity in the mixing channel.

. In addition, the corresponding average cell Pe numbers (Pe_Δ = uΔx/D) for different cases were also tabulated since this parameter helps us understand and control the amount of false diffusion in the system. Initially, numerical diffusion analysis was conducted using the smallest molecular diffusion constant in Table 2 (i.e., D_M = 3 × 10⁻¹⁰ m²/s) for all mesh levels of micromixer configurations (M1, M2, M3, M4, M5, and M6) and for two flow scenarios (Re = 240 and 120). Later, results obtained were evaluated and other molecular diffusion constants were tested for only mesh levels M1, M3, and M6 since false diffusion magnitude in transport simulations may be well defined based on the cell Peclet number, which carries the properties of grid size, velocity, and molecular diffusivity. Therefore, characterization of the numerical errors for three different grid sizes (or Pe_Δ) will provide sufficient information on the evolution of the false diffusion trend in the system.

5. Results and Discussion

5.1. Numerical Diffusion in Fluid Flow and Scalar Transport (D_M = 3 × 10⁻¹⁰ m²/s)

In numerical micromixer investigations, mesh characterization is a critical step to identify and control numerical errors which fundamentally originate from the properties of the mesh. In the current passive micromixer literature, there is no standard and accepted procedure for grid sensitivity analysis. Generally, tested mesh properties (e.g., grid size, grid type, mesh density, and element positioning etc.) are chosen as a rule of thumb or they are based on computational limits. However, these properties need to be considered and used appropriately depending on the physical problem type to avoid reporting unphysical mixing results in numerical simulations as remarked in References [9,10,11] and in also the present study. For instance, the selected mesh level in References [19,25] was primarily considered to resolve the flow field instead of solute transport although quite high Peclet numbers (e.g., on the order 10⁴–10⁶) were examined. In several papers [14,17,27,28], very close mesh densities were employed for grid studies which implies that the effects of numerical diffusion cannot be revealed clearly. In these applications, numerical simulations may be considered as grid-independent with a minimal error percentage even though a serious amount of error remains in the reported results. On the other hand, in some studies [2,13,14,17,27], low or moderate Re or Pe cases were used to determine a grid level, but higher Re or Pe case results were reported in these studies which are inconsistent.

Considering the above practices, six different mesh levels were determined for the micromixers designed to observe numerical diffusion in fluid flow and scalar transport. The density difference between M1 and M6 levels is around 14 million and 13 million for four-inlet and two-inlet micromixers respectively. Such a big density difference between meshes is necessary to capture the numerical diffusion effects in terms of the characteristics of the physical problem. Figure 5 shows the mesh study results for fluid flow at Re = 240 based on two different flow parameters, pressure drop in micromixers and velocity distribution at the exit of the mixing box (Line-I in Figure 3) where the most complex flow profile is observed. It is clear from Figure 5c,d that fluid flow was resolved identically for all mesh levels. The maximum difference between M1 and M6 for both parameters is less than 1%. As mentioned before, this small difference between grid levels occurred as a result of the relatively high kinematic viscosity (e.g., ν = 10⁻⁶ m²/s) of the fluids which leads to very low average cell Reynolds numbers (Re_Δ = uΔx/ν) in the computational domain. At Re = 240, even the biggest grid size (M6) results in an average Re_Δ number around 6 in the mixing channel which creates an insignificant amount of numerical diffusion in the flow solution. Meanwhile, the numerical diffusion term in fluid flow is used for numerical viscosity (or artificial viscosity) as discussed in Reference [11].

In contrast to the consistency between mesh levels in terms of the resolution of the flow field, transport simulations showed high discrepancy as a result of higher Pe_Δ numbers as given in Table 2, when D_M = 3 × 10⁻¹⁰ m²/s, which indicates the occurrence of sharp concentration gradients in the computational domain. Figure 6 and Figure 7 show scalar transport results for different mesh levels at Re = 240 and 120 flow conditions respectively. In each figure, while “a and b” and “c and d” plots represent four-inlet and two-inlet micromixer configurations, “a and c” and “b and d” plots show concentration distributions on Line-I and Line-II, as positioned in Figure 3, respectively. In both flow scenarios, while all mesh levels in both micromixer configurations exhibit relatively similar concentration distributions at the exit of the mixing box, these concentration trends differentiate at z = 500 µm in the mixing channel. This is because swirl motion starts in the mixing box and continually develops in the streamwise direction. Therefore, during the rotational flow of fluid pairs in the mixing channel, the transport solution starts producing numerical diffusion depending on the swirl profile that develops and the magnitude of the average Pe_Δ number for a specific mesh level. As a result of the flow patterns generated, concentration trends on the same sampling lines and variations between mesh levels are quite different for four-inlet and two-inlet configurations. In both micromixer configurations, however, there is a distinct difference between mesh levels M1 and M6 in terms of the resolution of the concentration field at z = 500 µm as shown in plots (b) and (d) of Figure 6 and Figure 7. Such a large variation emerged as a result of a doubled average Pe_Δ number between mesh levels M1 and M6 which are 10,000 and 20,000 for Re = 240 and 5000 and 10,000 for Re = 120 respectively. Meanwhile, variations between all mesh levels are obviously smaller in the Re = 120 when compared to the Re = 240 case, due to smaller average Pe_Δ numbers which generate relatively less numerical diffusion. On the other hand, although increasing the mesh density helps to resolve differences in concentration trends, still the finest mesh may contain a substantial amount of numerical diffusion because the average Pe_Δ number is still in the order of 5000 even for the best-case scenario tested in this section (e.g., Re = 120, M1 Level, and D_M = 3 × 10⁻¹⁰ m²/s).

To further investigate the numerical diffusion development in the micromixers, mixing on different cross-sections between the entrance and outlet of the mixing channel were measured and graphed as shown in Figure 8 in which “a and b” and “c and d” plots show four-inlet and two-inlet micromixer configurations and “a and c” and “b and d” plots show Re = 240 and 120 scenarios respectively. In the Re = 120 case, while mesh levels predict a similar amount of mixing values at the entrance of the mixing channel, measured mixing values are different due to developing mixing in the mixing channel. These differences between mesh levels emerge as a result of numerical diffusion during the mixing process in the mixing channel since each mesh level resolves different concentration profiles as previously shown in Figure 6 and Figure 7. Depending on the flow profile created and the magnitude of the swirls in the mixing channel, variation between mesh densities increases until a certain distance in the mixing channel is reached. After this point, mesh levels exhibit a mild convergent tendency. At Re = 240 scenario, however, while the two-inlet design shows a similar mixing estimation trend as observed in the Re = 120 case, the four-inlet configuration draws quite a different profile. As the variation between mixing indexes of different mesh densities increases until the z = 500 µm in the mixing channel, after this point, the variation declines and convergence is observed at the outlet of the micromixer. Thus, it should be noted that selecting the concentration sampling points for grid studies becomes important to reveal the actual contradiction between mesh levels. Using Equation (7), Figure 9 shows the comparison of mesh levels with the finest mesh in terms of mixing index on the x-y plane at different z-distances of the mixing channel. Similarly, while “a and b” and “c and d” plots show the results for four-inlet and two-inlet micromixer designs, “a and c” and “b and d” plots represent Re = 240 and 120 scenarios respectively.

$$\% Differenc{e}_{M(x)-M1}^i=\left|\frac{M{I}_{M(x)}^i-M{I}_{M1}^i}{M{I}_{M1}^i}\right|\times 100\hspace{1em}\{\begin{array}{l}i; z=500,1000,2000 \mathsf{\mu }\mathrm{m}\\ x; 2,3,4,5,6\end{array}$$

(7)

It is noticeable from Figure 9a–d that mesh refinement significantly reduces the numerical diffusion errors in both micromixer types and flow cases. Besides, the maximum difference occurs between M6 and M1 meshes as expected and gradually diminishes with increasing mesh density until the level of M2 and M1 is reached, which is around 10% for all cases. On the other hand, the mixing outcome obtained from separate locations results in different trends for mesh comparisons. Namely, discrepancies between meshes with the finest level reach the maximum at the z = 500 µm sampling point and beyond this point it starts decreasing across the mixing channel. While the two-inlet design reacts to mesh refinement and sampling regions similarly for both Re = 240 and 120 flow scenarios as shown in Figure 9c,d, four-inlet configuration shows quite low differences beyond the z = 500 µm sampling point at Re = 240 as shown in Figure 9a. Therefore, for this case, the use of the outlet mixing index in the mesh study may seriously mislead the evaluation of numerical diffusion because, as the difference between M6 and M1 is around 7% at the outlet, this is 42% at the z = 500 µm point. This contradiction is also observed in all other cases at relatively lower magnitudes. The discrepancy observed at different points and the converging tendency of different mesh resolutions across the mixing channel require an explanation. Figure 10a,b shows the fluid flow and scalar transport in the mixing channel of the four-inlet design for both Re = 240 and 120 scenarios respectively. As can be seen in Figure 10a, the swirl motion starts at the entrance of the mixing channel and continues strongly until z = 500 µm. After this point, the intensity of the generated swirl is dampened and fluid pairs flow along the mixing channel with a relatively smoother rotational movement. Hence, the maximum amount of numerical diffusion is produced between z = 100 µm and 500 µm depending on the grid resolution used and after this point the numerical errors generated are averaged and transported in the mixing channel as mixing.

The same explanation is also true for Re = 120 scenario, but the difference between results at z = 500 µm and beyond is less than that of Re = 240 case. This is mainly as a result of a relatively lower flow velocity and moderate swirl profile generated in the mixing channel as shown in Figure 10b. In the case of the two-inlet design, while the flow profile does not exhibit an even swirl form at Re = 240, the swirl profile generated is similar to the four-inlet design at Re = 120 as shown in Figure 11a,b respectively. Nevertheless, the maximum difference between mesh resolutions and the finest mesh is still observed at z = 500 µm for both flow scenarios. As discussed in Section 3, scalar transport with two-inlet injection is expected to produce less numerical diffusion since this configuration will create a relatively lower rotational contact surfaces between fluids. However, when four-inlet and two-inlet design outcomes are compared qualitatively at Re = 120, this is not as projected in the 2-D test case previously. Although fluids in the four-inlet micromixer create more contact surfaces at Re = 120 as shown in Figure 10b, the green regions in this figure are lower than that of the two-inlet solution as displayed in Figure 11b. If both figures are compared with respect to cross-sections at the entrance of the mixing channel (z = 100 µm), the two-inlet design shows much more green regions as opposed to four-inlet’s distinct blue and red color pattern. This is because the uniform velocity magnitude applied from each inlet of the two-inlet micromixer is two times higher than that of the four-inlet in order to provide Re = 240 and 120 flow conditions. Therefore, the flow regimes generated in the mixing boxes of both micromixer configurations are different as displayed in Figure 12a,b, which shows the central plane of the mixing box at z = 50 µm for both micromixer types and flow cases. It is clear from the figure that in contrast to the balanced four-sided fluid injection structure and lower inlet velocity of the four-inlet design, higher flow rates and two-sided inlet orientation generate a strong vortex inside the mixing box of the two-inlet micromixer for both flow conditions.

In addition, the strong vortex inside the mixing box creates two different swirls at the entrance of the mixing channel as one is at the center of the channel and the other one is around this central swirl as shown on different cross-sections of Figure 11a. It is obvious that while the green areas, which show a fully mixed state, on cross-sections at z = 100 µm are generated at the beginning of the mixing channel of the four-inlet micromixer, these green regions appear in-between two swirls and are carried from the mixing box of the two-inlet design as shown in Figure 10 and Figure 11 respectively. This difference may also be seen when the images in Figure 10, Figure 11 and Figure 12 are compared. As a result of higher inlet velocities in the two-inlet design, fluid bodies coming from both inlets encapsulate each other several times by rotating around the center of the mixing box as shown in Figure 13a,b. Besides, the size of the vortex profile created inside the mixing box exceeds the dimensions of the mixing box exit, at which the finest mesh elements are used as displayed in Figure 4. Hence, a higher average Pe_Δ number around the exit section and repeated mesh-flow disorientation inside the large vortex cause a drastic increase of numerical diffusion generation in the mixing box. In view of these results, false diffusion production is mainly controlled by the mixing channel and mixing box in four-inlet and two-inlet micromixer designs respectively. Accordingly, while the grid size distribution used inside the mixing box is a good strategy for a four-inlet micromixer, in the case of the two-inlet configuration, smaller mesh elements need to be positioned across the mixing box in order to control the amount of numerical diffusion produced.

5.2. Quantification and Analysis of Numerical Diffusion

In swirl-induced passive micromixers, higher flow rate requirement to create a swirling motion and very low diffusion constants inevitably lead to high Pe numbers. As it is shown in this study, the numerical solution of high Pe transport systems is quite challenging in terms of controlling the production of numerical diffusion throughout the system. Although the total mesh element numbers used are around 16.4 and 10 million for mesh levels M1 and M2 respectively, still, the discrepancy between these mesh levels is around 10% in terms of measured mixing. The effect of change in the numerical diffusion magnitude beyond M1 density is still not known due to the high computational cost. Considering the complexity of numerical diffusion analysis, Liu [10] developed a method to quantify the average numerical diffusion in steady flow systems utilizing the scalar gradient in numerical solutions and found that when average Pe_Δ = 50 is used, the numerical solution generated a negligible amount of numerical diffusion in the tested microchannel mixer. This method was later used in References [9,11] to analyze false diffusion in different micromixer types. In this technique, the combined effects of molecular diffusion (D_M) and numerical diffusion (D_N) are defined as effective diffusivity (D_Effective) as formulated in Equation (8). Therefore, Equation (8) makes it possible to evaluate numerical diffusion by setting the molecular diffusion constant equal to zero. In this case, effective diffusivity gives the average false diffusion amount in the numerical solution. On the other hand, if numerical diffusivity is zero in a system, effective diffusivity will be equal to the prescribed molecular diffusivity in the system. Effective diffusivity is calculated using a set of equations which may be found in the referenced articles in detail.

$$D_{Effective}\approx D_M+D_N$$

(8)

Table 3. Computed numerical (D_N, m²/s) and effective (D_Effective, m²/s) diffusion constants for different micromixer designs, flow scenarios, and mesh densities.

Mesh Level	Four-inlet Design				Two-inlet Design
	Re = 240		Re = 120		Re = 240		Re = 120
	DN	DEffective	DN	DEffective	DN	DEffective	DN	DEffective
M1	1.80 × 10−8	1.93 × 10−8	2.36 × 10−9	2.93 × 10−9	2.17 × 10−8	2.24 × 10−8	3.74 × 10−9	4.41 × 10−9
M2	2.62 × 10-8	2.74 × 10−8	3.64 × 10−9	4.30 × 10−9	3.11 × 10−8	3.19 × 10−8	5.44 × 10−9	6.17 × 10−9
M3	3.48 × 10−8	3.60 × 10−8	5.16 × 10−9	5.94 × 10−9	4.11 × 10−8	4.18 × 10−8	7.32 × 10−9	8.09 × 10−9
M4	4.41 × 10−8	4.52 × 10−8	7.13 × 10−9	8.04 × 10−9	5.21 × 10−8	5.27 × 10−8	9.40 × 10−9	1.02 × 10−8
M5	5.60 × 10−8	5.70 × 10−8	9.70 × 10−9	1.07 × 10−8	6.62 × 10−8	6.67 × 10−8	1.22 × 10−8	1.29 × 10−8
M6	6.59 × 10−8	6.69 × 10−8	1.24 × 10−8	1.36 × 10−8	7.83 × 10−8	7.89 × 10−8	1.47 × 10−8	1.54 × 10−8

shows the computed numerical diffusivity (when D_M = 0) and effective diffusivity (when D_M = 3 × 10⁻¹⁰ m²/s) constants from the numerical simulations of all micromixer designs, flow conditions, and mesh densities. First, it should be noted that numerical diffusion constants obtained show the unphysical diffusivity in the system for the designed micromixer types, flow conditions, mesh levels, and numerical scheme employed. Notably, estimated numerical and effective diffusivity constants are very close to each other for all scenarios, which indicates that the effect of molecular diffusion is severely masked by false diffusion. In both micromixer designs, effective diffusion constants calculated are one and two orders of magnitude higher than given molecular diffusion constant for Re = 120 and 240 flow cases respectively. Besides, increasing the mesh density also reduced the numerical diffusion generated for all tested scenarios. The two-inlet design simulations showed higher numerical diffusivities than that of the four-inlet micromixer structure especially in the Re = 240 flow scenario. This is obviously due to the strong vortex formation inside the mixing box as explained in the previous section. Nevertheless, all flow and mesh scenarios tested for both micromixer designs failed to recover the given molecular diffusion constants and exposed high numerical diffusion errors.

For a better presentation of Table 3, effective diffusivities were normalized by a molecular diffusion constant (D_M = 3 × 10⁻¹⁰ m²/s) and numerical diffusivities for each mesh level as shown in Figure 14a,b, respectively, in which the numbers above the trendlines represent the average Pe_Δ numbers for different mesh levels. The ratio of D_Effective/D_M describes the performance of the scalar transport simulation in terms of false diffusion production in the solution. Namely, if the amount of numerical diffusion generated approaches to zero, the molecular diffusivity used in simulations will be recovered by the computed effective diffusivity and, therefore, the ratio will be approaching to 1. Figure 14a clearly shows that the false diffusion amount produced in the Re = 120 scenario is significantly lower than that of Re = 240 for both micromixer types.

The two-inlet and four-inlet designs show a variation at same average Pe_Δ numbers which implies that the amount of mesh-flow misalignment inside the computational domain—hence the numerical diffusion generation tendency—is quite different in both micromixers. Accordingly, it should be noted that considering the cell Pe number alone in the control of false diffusion may yield misevaluation of the errors generated because it is obvious from Figure 14a that the two-inlet design is prone to create more numerical diffusion than the four-inlet design in the Re = 240 flow scenario. This difference increases when coarser grid elements are used in simulations. Even for the best-case scenarios tested in this study, the predicted effective diffusivities are approximately 10 and 15 times higher than the simulated molecular diffusivity for four-inlet and two-inlet designs respectively. As mentioned earlier, the smallest average Pe_Δ number is on the order of 5000 which is a large number to control false diffusion in numerical simulations.

Similarly, D_Effective/D_N may also be used for the false diffusion evaluation as shown in Figure 14b. In this case, when the ratio is one, the estimated effective diffusivity completely reflects the produced numerical diffusion in the simulation. When the produced false diffusion in a solution is very small, the ratio will be several orders higher than one. It is clear from the Figure 14b that even for the minimum average Pe_Δ number, the ratio is on the order of 1.25 which indicates that effective diffusivity is still close to the numerical diffusivity in the solution. When the two plots are compared, there is a consistency between computed effective diffusivities and numerical diffusion constants, which indicates that the employed method in Reference [10] coherently characterizes and quantifies the average numerical diffusion generation in a scalar transport simulation.

According to the outcome given in Table 3, the tested scenarios cannot reflect the physical effects of employed molecular diffusion constants as a result of quite high average Pe_Δ numbers, thus mixing values in Figure 8 are still masked by the numerical diffusion. Therefore, to observe the effects of smaller average Pe_Δ numbers, higher molecular diffusion constants were tested for mesh levels M1, M3, and M6, keeping the same flow conditions. Initially, the original molecular diffusion constants used (i.e., D_M = 3 × 10⁻¹⁰ m²/s) were increased 10 times for both flow conditions by which average Pe_Δ numbers were reduced 10 times. Later, considering the false diffusion production tendency of flow scenarios, 200 and 50 times higher molecular diffusivities were tested for Re = 240 and 120 cases respectively. The molecular diffusion constants used and the corresponding average Pe_Δ numbers may be found in Table 2. As displayed in Figure 15, which shows the change of the D_Effective/D_M ratio with respect to three different mesh levels and molecular diffusion constants, using a smaller average Pe_Δ number increases the recovery of the simulated molecular diffusion constant for both micromixer types and flow conditions. For instance, in the Re = 240 scenario, reducing the average Pe_Δ number from 10,000 to 1000 by using the D_M = 3 × 10⁻⁹ m²/s, decreased the ratio from 64 to 8 for the four-inlet design and 75 to 9 for the two-inlet design at mesh level M1 as shown in Figure 15a,c respectively. However, the observed ratios 8 and 9 at the M1 level are still very high to reflect the physical effects of the diffusion constant employed. The numerical diffusion generated at this level is on the order of 10⁻⁸ for both micromixer types as given in Table 3, which implies that numerical diffusion is one order of magnitude higher than the molecular diffusion constant and still dominates the physical diffusivity in the numerical solution. When D_M = 6 × 10⁻⁸ m²/s for the same flow scenario, the ratio reduces to 1.31 and 1.36 points at the M1 mesh level (Pe_Δ = 50) for the four-inlet and two-inlet micromixers respectively. In this case, the negative effects of numerical diffusion are mostly suppressed due to a tolerable average Pe_Δ number and, therefore, the ratio gets closer to one, which means that physical diffusivity is mostly reflected in simulations. For other mesh levels M3 (Pe_Δ = 70) and M6 (Pe_Δ = 100), the D_Effective/D_M ratio is “1.61 and 2.13” and “1.69 and 2.32” for the four-inlet and two-inlet micromixer designs respectively. Thus, when the average Pe_Δ is 100, the recovered diffusivity from the numerical solution is more than two times higher than the molecular diffusion constant used in the Re = 240 flow scenario. When Re = 120, the ratios at M1 level (Pe_Δ = 500, D_M = 3 × 10⁻⁹ m²/s) are 2.18 and 2.80 for the four-inlet and two-inlet configurations as shown in Figure 15b,d respectively. In addition, in the case of D_M = 1.5 × 10⁻⁸ m²/s, which results in an average cell Peclet number 100 at M1 grid level, the ratios obtained for the four-inlet and two-inlet micromixers are 1.34 and 1.42 respectively. If these ratios are compared with that of the Re = 240 scenario at the same mesh density M1, it is obvious that the physical diffusion recovery performance of both flow scenarios is almost equal while the Pe_Δ is two times higher in the Re = 120 scenario. Thus, when Re = 120, numerical simulations can tolerate higher average grid Peclet numbers in terms of numerical diffusion production, because swirling flow at Re = 120 is typically less when complicated to that of Re = 240 for both micromixer types. Consequently, while numerical diffusion generation is mostly determined by the magnitude of the cell Peclet number, the flow pattern created is similarly important since mesh-flow alignment of fluid pairs is controlled by the flow pattern.

In Reference [10], it is advised that when mixing is completed at very early stages in the micromixer, including the entire micromixer domain in effective diffusivity computation may result a substantially averaged and imprecise value. In this research, the maximum mixing index is observed at the outlet of the four-inlet micromixer in the Re = 240 flow scenario when D_M = 6 × 10⁻⁸ m²/s is used. The evolution of the mixing index values on x-y planes at z = 500, 1000, and 2000 µm in the mixing channel are 73, 95, and 98% respectively for the above case at the M1 level. Therefore, the effective diffusivity values computed mostly reflect the actual diffusivities in the computational domain since mixing is still developing through the outlet channel. In addition, concentration distributions at z = 500 µm (Line-II in Figure 3) for mesh levels M1, M3, and M6 are shown in Figure 16a,b for molecular diffusivities 3 × 10⁻¹⁰ and 6 × 10⁻⁸ m²/s respectively. As shown in Figure 16b, indeed, variations between mesh levels are significantly decreased because the above mesh densities resolved the transport domain at Pe_Δ numbers 50, 70, and 100 respectively. This figure visually reflects the computed D_Effective/D_M ratios. For example, when the ratios are 64 and 223 for mesh densities M1 and M6, respectively, the observed discrepancy is high between these two levels as plotted in Figure 16a. However, when the ratios drop to 1.36 and 2.32 by means of a higher molecular diffusion constant, Figure 16b shows that similar concentration trends are obtained for M1 and M6 mesh densities.

In this research, it was shown that characterization and quantification of false diffusion errors in swirl-induced mixing systems are critical to analyze the mixing performance and report physically reliable results. It should be noted that although numerical diffusion effects on physical mixing were analyzed in swirl-based mixing systems, the findings in this research are also valid for high Pe scalar transport systems in which secondary flows are dominant. In grid-based numerical techniques, diminishing the numerical diffusion errors to negligible levels at high Pe scenarios is possible when mesh density is increased. However, in specific cases, such as swirling flows, where grid-flow orientation is continuously violated, this approach may become unfeasible due to an extremely increased computational cost. For example, to render an average Pe_Δ number 10 in the mixing channel for the Re = 240 and D_M = 3 × 10⁻¹⁰ m²/s scenario, approximately 1 × 10¹⁵ hexahedron elements need to be used in the computational domain. Such a mesh density is obviously far beyond the computational capacity of today’s workstations. In this case, alternative approaches to the grid-based methods should be considered. In the current literature, several researchers proposed and practiced specialized particle-based numerical techniques for high Peclet transport systems, as can be seen in References [29,30,31,32,33,34] and the references therein. The algorithms developed in these studies present different computational advances which are beyond the scope of this study. However, these methods become prominent for elimination of the adverse effects of numerical diffusion errors at high Pe cases and less computational power requirement when compared to the conventional grid-based numerical methods.

6. Conclusions

In this study, 3-D swirl-induced micro-scale mixing systems were numerically investigated using an FVM-based CFD solver. Numerical simulations were conducted employing two different passive micromixers with two- and four-inlet fluid injection configurations. Fluid flow and scalar transport in these systems were studied in terms of numerical diffusion effects on the physical problem under different flow conditions, mesh densities, and molecular diffusivities. Before the 3-D micromixer tests, numerical diffusion effects for various scenarios were described using a 2-D test case. It was found that when the flow is exactly orthogonal to the grid boundaries, the scalar transport solution does not produce significant false diffusion. In contrast, oblique mesh and flow orientation cause unphysical diffusion in the numerical solution, the amount of which changes depending on the angle between mesh elements and flow. In addition, grid refinement and second-order discretization significantly reduced the produced false diffusion in the simulations.

In 3-D swirling micromixers, the flow domain at Re = 240 was similarly resolved by six different tested mesh densities with a maximum difference less than 1% between M1 and M6 levels as a result of very small average cell Re numbers (e.g., Re_Δ ≈ 6 for M6). While the four-inlet design generated a smooth swirl profile in the mixing box as a result of a balanced fluid injection structure, two-sided fluid injection and higher inlet velocities caused a strong vortex formation in the mixing box of the two-inlet design.

In the case of the scalar transport solution at Re = 240, the difference between mesh levels M1 and M6 increased up to 60% at z = 500 µm in the mixing channel because of a very large Pe_Δ which is on the order of 1 × 10⁴ when the tested molecular diffusivity is 3 × 10⁻¹⁰ m²/s. For the same molecular diffusivity, the measured difference between two fine meshes, M1 and M2, was still around 10% at both Re = 240 and 120. Thus, the false diffusivity amount in the numerical solutions were quantified for the characterization of numerical diffusion limits in both micromixer types. It was found that at Re = 240, recovered effective diffusivities from numerical solutions were 64 and 75 times higher than the tested molecular diffusion constant even at the M1 mesh level of four-inlet and two-inlet micromixers respectively. In the Re = 120 flow case, the above numbers dropped to 10 and 15 respectively. In all scenarios, the generated false diffusion amount in the two-inlet micromixer type was found to be higher than that of the four-inlet configuration as a result of a complex vortex profile in the mixing channel. It should be noted that for both micromixer types, a tolerable amount of numerical diffusion was generated when the molecular diffusion constants used were higher than the false diffusivity produced in the numerical solution. Namely, when Pe_Δ were 50 and 100 for Re = 240 and 120 flow scenarios respectively, the ratio of D_Effective/D_M was found to be around 1.35 which indicates that the recovered effective diffusivity is close to the tested molecular diffusion constant. Beyond these Pe_Δ numbers, physical molecular diffusivity effects in numerical solutions were masked by false diffusion and, therefore, the mixing values obtained show unphysical mixing in micromixers.

Consequently, it is shown that numerical diffusion generation in 3-D swirling passive micromixers depends on cell Peclet magnitude and the flow pattern formed. When cell Pe number is larger than a certain value and grid-flow alignment is constantly disturbed, numerical schemes cannot resolve high scalar gradients accurately and produce unphysical numerical diffusion in the solution. Thus, numerical simulations of advection dominant swirling systems need to be conducted carefully and the numerical diffusivities generated should be quantified to avoid reporting unphysical mixing results.