Design and validation of microfluidic parameters of a microfluidic chip using fluid dynamics

The internal fluidic parameters of microfluidic channels must be analyzed to solve fundamental microfluidic problems, including microscale transport problems involving thermal analysis, chemical reactivity, velocity, pressure drop, etc., for developing good-quality chemical and biological products. Therefore, the characterization and optimization of the interaction of chemical and biological solutions through microfluidic channels are vital for fluid flow design and engineering for quality assurance in microfluidic platforms. As the internal structures and kinetics of microfluidic channels are becoming increasingly complex, experiments involving optimal fluidic and transport designs are challenging to perform with high accuracy. However, highly integrated simulation tools can guide researchers without specialized computational fluid backgrounds to design numerical prototypes of highly integrated devices. In this study, a microfluidic chip with two inlet wells and one outlet well was fabricated from polydimethylsiloxane following which simulations were performed using an ANSYS Fluent tool influenced by computational fluid dynamics at a nearly identical scale. The pressure drop and velocity profiles of the interaction of two pH buffer solutions (pH 4 and 10) through the designed microfluidic chip were qualitatively estimated from experimental data analysis and validated with the simulation results obtained from the CFD-influenced ANSYS Fluent tool. © 2021 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0056597


I. INTRODUCTION
The testing of chemical and biological samples through microsized devices has attracted special attention for reliable, early, and quick medical diagnosis and quality assurance in microfluidic platforms. Therefore, the intensive fluid dynamics-based analysis of liquid samples through such micro-devices that provide efficient output can ensure substantial health benefits for the next generation. Transport phenomena 1,2 on the microscale have received special attention owing to increasing interest in efficient and sustainable processes. Despite the laminar flow conditions in the interaction of different fluids on the microscale, numerical simulations with high accuracy are still necessary, in addition to the experimental visualization of flow of mixing fluids according to the simulation scale, which is a basic necessity. In reality, the exact data of geometries and wall conditions of internal microflow channels and data on chemical media, such as diffusion coefficients and reaction rates, are not published before performing the simulated and experimental visualization of the flow fields of fluids. Therefore, to understand and optimize microflow channels as well as to trace the pattern of internal microfluidic parameters, microscale flow must be visualized in the simulation scale first and then verified through experimental visualization according to the simulation scale. Since the early 2000s, researchers in fluid mechanics have relentlessly visualized the flow fields of fluids to trace the pattern of internal fluid flow paths with the objective of achieving a high spatial resolution using moderate optical arrangement methods. 1, 3 Ahmed et al. 4 proposed, analyzed, and visualized the microscale flow field of the interactions of two different pH buffer fluids by using terahertz (THz) image sensing technology to trace the flow patterns with respect to the change in the intensities of peak values of THz amplitude data. ANSYS Fluent software influenced by computational fluid dynamics (CFD) 5-7 is a particularly effective tool for the study and analysis of microfluidic flow behaviors. This computational analysis method has many merits, including robust device design and the ability to simulate coupled and complex physics quickly at a low cost. 8,9 When two or more different viscous fluids mix through microflow channels with variations in velocities, it is necessary to visualize and analyze the flow patterns of fluids at each point of the internal flow channels by characterizing the nature of fluid flow and the different patterns of internal dimensions of the channels. Therefore, we need to conduct 2D simulations using the ANSYS Fluent simulator to obtain simple and precise data reflecting fluid-flow patterns under the influence of friction due to variations in the speeds of different fluids having different viscosities within the internal layers of microfluidic channels. Thus, performance evaluation and optimization with different dimensions of internal microfluidic channels can be achieved by analyzing the characteristics of internal microfluidic parameters, such as pressure drops, velocities, densities, and viscosities.
The digital pressure machine (DPM)-controlled fluid flow through the capillary design concept-based microfluidic chip is schematically shown in Fig. 1. In the present study, pH 4 and 10 buffer solutions were injected through two inlet wells of a newly designed microfluidic chip by using a digital pressure machine (Model fusion 710, Chemyx Inc.), and the mixture of these two solutions was deposited into a beaker through the only outlet well of the chip, as shown in Fig. 2. From the deposited mixture of the two pH buffer solutions, the fluid flow rates were calculated volumetrically from the total measurement time. As part of the performance evaluation of the chip, pressure drops and velocity profiles from the digital pressure machine were estimated at the two inlet wells and the outlet well. Subsequently, from the experimentally measured fluid flow rates at the outlet well, curves of velocity vs fluid flow rate and pressure drop vs fluid flow rate were drawn to illustrate the continuity equation and Hagen-Poiseuille equation of fluid dynamics. Finally, the experimental results were validated with simulated data obtained from a CFD-influenced ANSYS Fluent tool.

II. METHODS
Ethics approval is not required.

A. Preliminaries
According to the electrical circuit model, 10 to deduce the fluid flow rate, the microfluidic resistance (MFR) values and the given pressure drops applied to a microfluidic chip are used, followed by the Hagen-Poiseuille equation 11,12 where R f I/O is the microfluidic resistance (MFR) through each of the inlets and the outlet of microfluidic channels and Q I/O is the volumetric fluid flow rate through each of the inlets and the outlet. ΔP I/O is the pressure drop difference across each of the inlets and the outlet. From solutions of Navier-Stoke's theorem, circular or cylindrical and hydrodynamic resistances were calculated (Appendix A).

B. Experiments
Polydimethylsiloxane (PDMS) 13 is widely used for fabricating micro total analysis systems 14,15 and lab-on-chips. 16,17 In this study, a Sylgard ® 184 silicone elastomer base agent and a Sylgard ® 184 silicone elastomer curing agent as an adhesive were employed, and these materials were manufactured by Dow Corning (Midland, MD, USA). From these base and curing agents, we developed a PDMSbased microfluidic structure through an irreversible assembly process. For this purpose, a 3D structure of a microfluidic device with two inlet wells and one outlet well was drafted using SolidWorks software [Figs. 3(a) and 3(b)], translated into the STL format for 3D printing, and printed using a high-definition 3D printer (Agilista-3100; serial no.: 96M14458) in our laboratory [ Fig. 4(a)]. The base and curing agents, at a mass ratio of 10:1, were cross-linked to manufacture the designed PDMS microfluidic chip. During the fabrication process, the cross-linked base and curing agents were poured into the 3D-printed structure with a cross section of 20 × 20 mm 2 and a height of 7 mm from the base surface to fabricate a replica mold of the 3D structure [ Fig. 4(a)]. To confirm whether the replica mold had sufficient hardness, the poured solutions were placed in an oven controlled at 40 ○ C for 24 h. The replica mold was manually cut from the 20 × 20 mm 2 structure in dimensions of 10 × 10 mm 2 to manufacture the PDMS microfluidic chip. 4 The heights of all inlets and the outlet are 7 mm from the base surface and that of the internal microflow channels is 4 mm from the base surface [ Fig. 3(b)]. We attached the bottom part of the microfluidic chip to a glass substrate to ensure the intensive flow of fluids through the microflow channels of the microfluidic chip. 4 A top plate was screwed and tilted mechanically with a bottom plate by sandwiching the microfluidic chip along with the glass substrate between both plates [ Fig. 4(b)].
The viscosities of the pH 4 and 10 buffer solutions were experimentally measured as 0.5 and 0.59 mPa s, respectively. The viscosity of the pH 4 and 10 mixed solution is 0.58 mPa s. Each of the two syringes connected to two tubes was 2 m in length and 2 mm in diameter. The mixed accumulated buffer solution was transferred into the common junction well 4 and flushed out through the outlet well into a beaker by using a tube 1 m in length and 2 mm in

ARTICLE
scitation.org/journal/adv diameter. The Reynolds number, 11,12 Re, was estimated at the outlet well by calculating the average velocity from the experimentally measured values of different fluid flow rates (Appendix B). During the experiment, the digital pressure machine was set at different fluid flow rates by pumping two syringe-connected tubes automatically to allow the flow of pH 4 and 10 buffer solutions through the microfluidic chip. Different pressure drops corresponding to five different fluid flow rates set using the digital pressure machine were calculated for the two inlet wells from Eq. (1) by substituting the cylindrical resistance values calculated using Eq. (A1) for the two inlet wells of the pH 4 and 10 buffer solutions. The cylindrical resistances for the pH 4 and 10 solutions were calculated as 97.4 and 114.9 MPaΔsΔm −3 , respectively, using (A1) and the different dimensions from Figs The average pressure drop in each interval of the fluid flow rate (Fig. 5) is higher than that for each case at the outlet well in Fig. 6. This is because the syringe-connected tubes experienced a huge pressure drop as the air did not fully escape from inside the syringes.  To ensure that the digital pressure machine automatically pushes the plungers of the syringes, the air inside the plungers of the two syringes was compressed, creating a much higher pressure inside the two syringes. Consequently, this high pressure drop was automatically applied to the tubes connected with the syringes, which affected the effective diameters of the sophisticated structure of the two tubes instantaneously to handle the high pressure drop inside the tubes. Therefore, the two syringe-connected tubes accessed the two inlet wells by creating a higher overall impact on the pressure drop inside the two inlet wells than on the pressure drop of the outlet well. Additionally, there is another vital cause for such a higher pressure drop across the two inlet wells. The designed diameter of each of the two cylindrical inlet wells (1.1 mm each) was lower than that of the outlet well (2 mm). Consequently, the pressure drop across the two inlet wells (Fig. 5) was higher than that of the outlet well ( Fig. 6) according to Eq. (1), where the pressure drop is proportional to the fluid flow rate in each measurement. The average velocities, pressure drops, and measured average fluid flow rates at the outlet well are listed in Table I for reference. As shown in Fig. 7, the average velocity at the outlet well was estimated statistically by repeatedly taking five measurements of  fluid flow rates at the outlet well with respect to each of the set values of the fluid flow rate at the two inlet wells. To perform this measurement using the trial and error method, mixed pH 4 and 10 buffer solutions were flushed out automatically into a beaker through a tube connected with only the outlet well of the microfluidic chip. Subsequently, the volume of the flushed-out mixed fluids was calculated from the total occupied volume of the mixed fluids within the beaker in milliliters for a certain estimated total measurement time (32 min for each measurement). Thus, the fluid flow rates were calculated and then averaged. From the five readings of each fluid flow rate, five velocity values were calculated for each case, and the average velocity was calculated using Eq. (A4). The plot of the average velocity at the outlet well against the average fluid flow rate, shown in Fig. 6, indicates that the average velocity increased linearly with the average fluid flow rate with the use of the trial and error method during the experimental measurement time. For the data analysis, Appendix B is referred.

A. Physical model and SolidWorks modeling
The geometric details of the structure described in Figs. 3(a) and 3(b) were used as input data to simulate the fluid flows. The internal fluid flow parameters and channels were numerically simulated using the academic version of the ANSYS Fluent workbench. 18 A SolidWorks model was drawn in a Navier-Stokes design modeler with a geometric shape similar to that of the fabricated microfluidic channels of the designed microfluidic chip. The model was reduced based on the axial symmetry of the domain into a planar twodimensional section to fix better convergence criteria and reduce the computational cost. In the numerical model, the viscosities and densities of the two pH buffer solutions were considered, with other properties set to those of water, as listed in Table II. All microfluidic parameters given as inputs were calculated from experimentally measured data.

B. Governing equations
The multiphase volume of fluid (VOF) 19,20 approach with two Eulerian phases was considered to track the interface between the two pH buffer solutions. Moreover, to activate the dynamic coalescence behavior of the two pH buffer solutions in a stagnant environment, the continuity, energy, and momentum governing equations 18,21 of fluid dynamics were adopted in this numerical simulation. The Euler form of the equations was used (Appendix C).

C. Numerical settings
Using the CFD package 22 of ANSYS Fluent, a powerful commercial tool, we performed a numerical simulation with a residual convergence condition of 0.001 for pressure, density, and momentum and of 10 −8 for energy. We used a pressure-based solver with the absolute velocity for the transient time-dependent model. The fluid flow was assumed to be laminar. The boundary walls were imposed with no-slip conditions and set to a constant temperature of 300 K. As inputs, the applied boundary conditions were the velocity in the two inlet wells for both the pH 4 and 10 buffer solutions, which had the same given velocity and the different densities, viscosities, etc., as listed in Table II. The fluidic properties of both pH buffer solutions were assigned to a new material section in the Fluent database. The mixed pH buffer solution was considered in the pressure outlet. A simple algorithm for pressure and velocity coupling was used in the solution methods. Moreover, to achieve a spatially discretized least-squares cell-based gradient with PRESTO, the pressure and second-order upwind momentum for a first-order implicit scheme were added to the solution methods. The under-relaxation factors in solution controls were taken as half of the default to improve the convergence rate. Hybrid initialization was performed with time steps of 0.001 s. In the computational approach, fixed-type time advancement was assigned with 20 time steps for a maximum of 100 iterations. However, the solution converged within 450-800 iterations when the residual reached the limiting point.

D. Mesh dependency and validation of simulation
The simulation was validated with the experimental results of pressure drops in the designed microflow channels of two inlet wells and one outlet well with errors ranging from ±35% to ±4%. Using element sizes of 4.5 × 10 −2 , 5 × 10 −2 , 5.5 × 10 −2 , 6 × 10 −2 , 7 × 10 −2 , and 8 × 10 −2 mm with ten inflation layers at the boundary, meshes with an orthogonal rectangular structure were created for the model. All the simulations were then conducted for a 5 × 10 −2 mm element size with 10 185 elements, which yielded the lowest error based on trials. Tables III and IV present data validation and estimation using ARTICLE scitation.org/journal/adv  linear regression, respectively. Figures 8 and 9 graphically show the model validation and mesh dependency, respectively.

IV. RESULTS AND DISCUSSION
When two pH buffer fluids in microfluidic channels were flushed out through the outlet well, the outlet well itself was opened  at the atmospheric pressure and outside temperature. With the velocity of the pH buffer solutions mixed via the two inlet wells, the net mixed product of these fluids was estimated through the open atmospheric pressure of the outlet well. Therefore, in reality, the final fluidic parameters at the outlet were based on the mixture of the two pH buffer solutions, and the mixture was influenced by the external atmospheric pressure and the three-dimensional shape of the outlet well, the diameter of which was higher than that of the two inlet wells. However, in the case of two-dimensional simulations, we only considered fluids flowing through a planar surface without accounting for the effect of fluids flowing through three-dimensional microfluidic channels. Moreover, the CFD-influenced tool was an ideal, external pressure free environment to perform the simulation. As a result, the error seemed to vary from ±4% to ±35% in the two-dimensional simulation. The error percentage of the change in pressure drop in between the experimental and simulation results can be calculated as follows: However, from the experimental measurement, for instance, at a total fluid flow rate of (0.1 + 0.1) = 0.2 ml/min through both the inlets (Table I), we obtained a corresponding net fluid flow rate of 0.17 ml/min at the outlet. This is a result of the loss of a significant amount of fluid inside the surface walls of the microflow channels due to the consideration of a certain measurement time by instantaneously stopping the flow of fluids. Therefore, it was necessary to apply the linear regression method by taking the average of each pressure drop at the outlet well to overcome such practical difficulties. Considering the lowest error during the experiment, the mesh dependency test was applied for 10 185 elements for velocities ranging from 0.26 to 0.38 mm/s, as this number of elements yielded the lowest error in the trial process. In this velocity range, the minimum error from the experimental results was considered to validate the simulation results. If we consider the validation based on the average pressure drop using the linear regression method, from Table III, Fig. 10, from inlet points via intermediate points C 1 and C 2 and the common junction point C 3 to the exit point E. According to the fundamentals of fluid dynamics, as two fluids do not mix by flowing through their own layers separately, the pH 4 and 10 buffer solutions showed laminar fluid flow. Two different viscous solutions traveled through the robust internal surface of the designed microflow channels, and friction 24 occurred at the internal boundary walls of the microfluidic channels, which in turn created a tendency of slight deflection in speed 4 from that of the continuous flow of the two fluid flows. The higher the viscosity of a fluid, the higher is the tendency of the speed of fluid flow to deflect with increasing fluid flow rate as observed. Therefore, the yellowish and red fluid flow profiles in Fig. 10 were for the pH 10 solution, the viscosity of which (0.58 mPa s) is higher than that of the pH 4 solution (mapped with the surface of its flow path).
For each VF, the fluid flow paths showed a slight change at different velocities to provide a mixture of fluids at a velocity of 0.096 mm/s at the points A and B. Similarly, the same input velocity of 0.064 mm/s, corresponding to a fluid flow rate of 0.1 ml/min for the two pH solutions, was applied to the points A and B for different VFs. However, in the simulated data, there was no appreciable change in the fluid flow paths for different VFs. This indicates that the velocity profiles remained nearly constant as the VF increased. Therefore, the simulated data with an optimum fluid flow rate of 0.1 ml/min in Fig. 10(a) would be effective to achieve nearly ideal fluid flow patterns through the microflow channels because this fluid flow rate yielded consistent fluid flow behavior with almost the same laminar velocity profiles despite the increase in VF, which was our ultimate goal. In Fig. 10(b), with the increase in velocity at the two inlets according to the increase in VF, the fluid flow path was notably affected within the CFD-influenced environment.

V. CONCLUSIONS
Microfluidics technology has been proven to be an effective and powerful tool to connect engineering fields with the life sciences. In this study, we designed a microfluidic chip with two inlet wells and one outlet well and studied the flow of two pH buffer solutions through it. We proposed a 2D simulation approach to validate the experimental results, quantify the deviations between experimental and simulation results, and analyze the practical causes of such deviations. However, the percentage of deviation was at an expected level, as external atmospheric factors outside the microfluidic channels affect the overall output during the experiment. Meanwhile, simulated data were obtained in a CFD-influenced environment, which was more compact than that of the experiment. In the analysis of fluid flow patterns inside the microflow channels, we observed that it is more convenient to simulate calmer fluids and analyze their internal fluidic parameters.
The above analysis based on the proposed microfluidic chip can reveal the traces of interactions of flowing fluids, such as chemical solutions and abnormal biological cell division inside tissue samples and blood samples. Therefore, the designed microfluidic chip can help to detect pH imbalances in the human body due to the malfunctioning of metabolism, respiration processes, and cancer cells, enabling early medical diagnosis. Moreover, this analysis is expected to be effective for the pH control of foods in food industries, tablet coatings in pharmaceutical industries, and so on for quality assurance; pH control is highly important in these industries.
In the future, we will investigate a 3D simulation approach using the ANSYS Fluent simulator and simulate microfluidic chips ARTICLE scitation.org/journal/adv with more complex structures using terahertz laser-scanned images of microfluidic internal reactions. Such work may reveal details of chemical reactions inside the flow channels of complex microfluidic structures to expand the applicability of microfluidics technology. Moreover, to make further contributions in the field of biomedical engineering, we plan to fabricate microfluidic chips through which more than two chemical solutions can flow by conducting fluid dynamics-based analysis. Comparing the simulated and experimental results for the microfluidic chip with a conventional height of 100 μm is under way as the next step.

APPENDIX A: MATH AND EQUATIONS
For a laminar, isothermal, incompressible, and isotropic Newtonian fluid flow through a tube or well with a circular cross section and no-slip boundary conditions on the sidewalls of the microfluidic structure, the microfluidic resistance can be modeled 11 mathematically as follows: where μ is the viscosity of the fluid, L is the length of the tube/well, and d is the diameter of the tube/well. 1 d 4 term indicates that the resistance values change vastly even with a small change in d. For a rectangular cross section, the fluid flow rate, Q R , can be formulated in terms of the pressure drop ΔPR over the microfluidic channel length L according to Eq. (1), The pressure drop, ΔPR, can be expressed as a function of the hydrodynamic resistance in the rectangular channel, R f H , which is also called the chip resistance, R f chip ; the fluid flow rate is similar to voltage drops in electrical current, where w 1 is the width of the rectangular channel and h is the height of the rectangular channel. Equation (A3) calculates the hydrodynamic or chip resistance. 12 The continuity equation 25 is as follows: where A I/O is the area of either the circular or rectangular path and V I/O is the average velocity of the fluid through each of the inlets and the outlet. The area of the inlets/outlet A I/O = 2πr(r + h), where r is the radius of the cylindrical-shaped inlets/outlet and h is the height. The Reynolds number Re is a dimensionless number that qualitatively characterizes a fluid in motion. It is expressed as the ratio of convective forces to viscous forces, where V I/O is the average fluid velocity in each of the inlets and outlets and D hh is the characteristic length of the micro-flow channel.

APPENDIX B: EXPERIMENTAL AND SIMULATION DATA CALCULATION
Re was estimated at the outlet well of the microfluidic chip by calculating the average velocity from the experimentally measured values of different fluid flow rates. In the outlet well, the density of the mixed solution was 0.9901 g/ml. According to the different experimentally measured fluid flow rates, the average velocity at the outlet well is calculated as where v is the fluid velocity (m/s), ρ is the fluid density (kg/m 3 ), p st is the static pressure (Pa), μ is the dynamic viscosity (mPa s), α 1 is the thermal diffusivity, T is the temperature, λ is the thermal conductivity, h is the internal energy, and Fcss is the surface tension. For the two phases, in each control volume, the density, viscosity, and thermal conductivity were calculated using the following equations: In general, for an n-phase system, all other properties in the volume-fraction-average method take the following form: x = ∑ αqxq. (C8) The VOF model treated the energy, Ee, and the temperature, T, as mass-averaged variables as follows:

ARTICLE scitation.org/journal/adv
where Eq for each phase is based on the specific heat of that phase and the shared temperature. In this model, the continuum surface force (CSF) model was used to calculate the surface tension with a non-conservative formulation.

DATA AVAILABILITY
The data that support the findings of this study are available within Appendix B of the article.