Accelerated Computational Fluid Dynamics Simulations of Microfluidic Devices by Exploiting Higher Levels of Abstraction

The design of microfluidic devices is a cumbersome and tedious process that can be significantly improved by simulation. Methods based on Computational Fluid Dynamics (CFD) are considered state-of-the-art, but require extensive compute time—oftentimes limiting the size of microfluidic devices that can be simulated. Simulation methods that abstract the underlying physics on a higher level generally provide results instantly, but the fidelity of these methods is usually worse. In this work, a simulation method that accelerates CFD simulations by exploiting simulation methods on higher levels of abstraction is proposed. Case studies confirm that the proposed method accelerates CFD simulations by multiple factors (often several orders of magnitude) while maintaining the fidelity of CFD simulations.


I. INTRODUCTION
A microfluidic device, or a Lab-on-a-Chip (LoC), is a device that performs lab operations on the microscale through a set of fluid manipulations [1].Such devices are commonly used for, e.g., personalized medical care [2], point-of-care diagnostics [3] (well-known examples of such devices are pregnancy tests [4] or the SARS-CoV-2 tests [5]) and the food industry [6].They have recently been proven to be more widely applicable in, e.g., geosciences [7] or fuel cell technology [8].In that regard, much potential lies in microfluidic devices.
The development of those microfluidic devices is a cumbersome and tedious process that often requires multiple expensive and time-consuming design cycles [9].To advance the design of microfluidic devices, reliable and quick simulation methods are necessary to predict whether a design works as intended.The simulation of Newtonian flow through a single microchannel can be performed using the Hagen-Poiseuille Law [10] and the flow profile across the channel can be assumed to be parabolic [11].However, channel-based microfluidic devices consist of multiple interconnected microchannels, which require more extensive simulation methods.
To this end, most designers utilize methods from Computational Fluid Dynamics (CFD), such as the Finite Volume Method (FVM, [12]), Finite Element Method (FEM, [13]), or the Lattice Boltzmann Method (LBM, [14]), to obtain results of good fidelity.In recent overviews for microfluidics modeling, the FVM [15], FEM [15], and LBM [16] are listed as numerical approaches to solving the Navier-Stokes Equations (NSE) for microfluidics.Another recent work describes the workflow of setting up simulations for microfluidic devices [17] with OpenFOAM v9.0 [18], which uses the FVM.Hence, these simulation methods can be considered state-of-the-art for modeling microfluidic devices.
However, CFD simulations can take up to days or even weeks, even on dedicated workstations [19].In practice, this obviously limits the use of CFD simulations for microfluidic devices to single components of the device [20].To fully understand and design a harmonious device, it is critical to simulate the behavior of all components and their interaction with each other.
Alternatively, high abstraction simulation methods, i.e., simulation methods that abstract the underlying physics on a higher level (also known as reduced-order modeling), generally simulate microfluidic devices almost instantly (i.e., in less than a second) [20].An example of such a method is to draw an analogy between the Hagen-Poiseuille law and Ohm's law and apply analogous methods from electrical circuit engineering to channel-based microfluidic devices [21].This approach is not limited to basic fluid flow, but can also be applied to problems with, e.g., droplets [22], or capillary-driven flows for paper-based microfluidics [23].Such methods may not provide results of comparable fidelity, but can still simulate some parts of microfluidic devices, e.g., channels, with relatively good accuracy [20].
In this work, the nature of these high abstraction simulation methods is exploited to substantially accelerate CFD simulations of microfluidic devices with hardly any loss of fidelity.To this end, a two-stage approach is proposed: First, regions of the microfluidic device are identified that can sufficiently be simulated at a high level of abstraction.Afterward, the corresponding simulations are conducted, and the respectively obtained results are communicated between the simulation engines.Case studies (using continuous channel-based microfluidic devices as a representative) confirm the potential of this approach: The proposed approach does not only generate simulation results faster than the original CFD method, but constantly does so by several factors or even several orders of magnitude-while, at the same time, maintaining the fidelity of the results.
The remainder of this work is structured as follows: First, we review the simulation methods for microfluidic devices in more detail, focusing on methods on low and high abstraction levels.Afterward, the accelerated CFD simulation method is described in a general fashion in Section III.Implementation details of this method are then provided in Section IV.Finally, we demonstrate the resulting solution and compare it against solutions of CFD simulations for a set of test cases in Section V-confirming that the accelerated method is faster than the CFD method for all test cases while maintaining the fidelity.Finally, the paper is concluded in Section VI.

II. SIMULATION METHODS FOR MICROFLUIDIC DEVICES
Simulation methods for microfluidic devices can be categorized into methods based on Computational Fluid Dynamics (CFD, [12]- [14]) on the low abstraction level, and methods on high abstraction levels [21]- [23], sometimes referred to as 1D methods.CFD methods can be considered state-of-the-art in the design of microfluidic devices [15]- [17], whereas the high abstraction level methods are often used to derive initial estimates during the design process [9].In this section, we review the CFD and high abstraction level simulation methods.

A. Review of CFD Methods
A low abstraction simulation of microfluidic devices can be obtained through methods from Computational Fluid Dynamics (CFD, [12]- [14]).In CFD, the fluid behavior is modeled by the Navier-Stokes Equations (NSE), which are the fundamental governing equations for fluid dynamics.We restrict ourselves to the incompressible NSE [12]- [14], which are given by the mass equation where u is the flow velocity vector, and the momentum equation where ρ is the fluid density, p is the pressure, and ν is the kinematic viscosity.
The analytical solution of the incompressible NSE exists only for a few simple fluid dynamics problems.For microfluidic devices, the equations generally have to be solved numerically.To this end, a wide variety of numerical methods has been developed.For example: • The Finite Volume Method (FVM, [12]).The FVM splits the computational domain into grid cells, and the NSE are solved numerically on each grid cell.By this, averaged values for the flow velocity and pressure are obtained in each cell.• The Finite Element Method (FEM, [13]).Similar to the FVM, the domain is split into cells.However, the solution is represented by a set of elements, e.g., polynomials.• The Lattice Boltzmann Method (LBM, [14]).Rather than solving the NSE, this method solves the Boltzmann equation-which can be proven, through the Chapman-Enskog theory, to macroscopically solve the NSE [14].CFD methods have in common that they, in a general sense, acquire results of good fidelity and are, therefore, helpful for the design of delicate microfluidic devices or components.
However, CFD methods require significant computational resources-in terms of computational memory and time.A CFD simulation of entire microfluidic devices can be hard to get right, and compute times can get up to days or weeks, even on dedicated workstations [19].

B. Review of Methods on High Abstraction Levels
Methods with a high level of abstraction simulate microfluidic devices by abstracting the underlying physics on a high level.Generally, these methods are based on solutions that can be obtained analytically for fairly simple problems.They map the results to similar problems under a set of simplifying assumptions, or models are obtained empirically through fitted data from experiments or pre-simulations [24].These methods are generally of poor fidelity, but results for entire microfluidic devices are usually acquired instantly (i.e., in less than a second) [20].More precisely, high abstraction simulation methods have been proposed for the following microfluidic platforms [25]: • Continuous channel-based microfluidics.This platform consists of a network of rectangular channels with width and height in the order of micrometers.Liquid flow through these channels can practically always be assumed to be laminar, and the flow profile in a channel can be accurately solved using the Hagen-Poiseuille law, i.e., where Q is the flow rate and R H is the hydraulic resistance of a channel.Using the Modified Nodal Analysis (MNA, [21]), the pressure and flow rates of all channels in a connected network can be calculated.• Droplet-based microfluidics.This platform has a network of channels, similar to the continuous channel-based microfluidics platform, with additional droplets of a fluid that is immiscible with the carrier fluid (continuous phase).The hydraulic resistance in Eq. ( 3) can be split into the resistance of the channel R channel H and the resistance of a droplet R droplet H present in that channel, i.e., Based on this and the MNA, this platform can be simulated on a high abstraction level [22].• Paper-based microfluidics.In paper-based microfluidics, a liquid is transported through a two-dimensional sheet of paper using capillary force.The one-dimensional transport of a liquid front through a porous medium is given by the Washburn equation, i.e., where L is the traversed distance of the fluid front, γ is the effective surface tension, D is the diffusivity coefficient, t is time, and µ is the dynamic viscosity.Based on Eq. ( 5), the capillary transportation of fluid can be simulated for porous channels with arbitrary cross-sectional shapes [23].

III. ACCELERATING CFD SIMULATIONS
Motivated by the fast computation time of methods with a high level of abstraction, the possibility to accelerate CFD simulations by exploiting said methods is investigated.In this work, we aim to accelerate CFD simulations for steady-state flow of continuous channel-based microfluidic devices as a representative.Firstly, the semantics of continuous channel-based microfluidics are covered in this section, and the potential for a faster simulation is highlighted.Afterward, a method will be proposed that explicitly exploits that potential and accelerates CFD simulations.Based on that, implementation details for this method and a summary of corresponding case studies, including evaluation results, are provided in the following sections.

A. Continuous Channel-Based Microfluidics
In continuous channel-based microfluidics, we consider a network as sketched in the middle of Figure 1.This example network has three inlets and one outlet and contains a homogeneous fluid (no mixture).It is depicted here as a two-dimensional network of channels with width w, and the extension to the third dimension for real-world microfluidic devices can be performed by adding a height parameter h.Furthermore, without loss of generality, we assume an adiabatic system and ignore gravity effects, such that the only relevant fields are the pressure p and velocity u of the fluid.However, as in this setup, the pressure and velocity fields of a fluid in a microfluidic network are generally complex and require dedicated methods to solve Eqs. ( 1) and ( 2).This is sketched in Fig. 1 for the pressure and velocity field of the steady-state flow at the location where channels cross (the red circles at the top).The contour lines of the pressure field and vectors of the velocity field are chaotic and not easy to predict.On the other hand, if we look at the pressure and velocity fields of a straight channel section (the green circles at the bottom), we notice a more organized and streamlined flow.The contour lines of the pressure field are straight, evenly spaced, and perpendicular to the channel.We can see here that, for straight channels, the pressure is a point-value along the channel and equal over the channel's cross-section.The velocity vectors are also organized and show a parabolic (for Newtonian fluids [20]) flow profile, as can be found, for a two-dimensional channel, by using where x is the direction parallel to the channel, and y is perpendicular to the channel.This is the cornerstone, on which Eq. ( 3) and, therefore, the high abstraction simulation method for continuous channel-based microfluidics is based.

B. The Proposed Accelerated Method
Based on the observation above, we can exploit the high abstraction level simulation method in regions where the flow is highly organized (green circles in Fig. 1) and use CFD simulations for regions where the flow is chaotic and hard to predict (red circles in Fig. 1).Exploiting the high-speed simulation feature of methods on a high abstraction level not only significantly reduces the compute time for large microfluidic devices (resulting in more favorable scaling of simulations) but also reduces the required memory.To this end, two steps must be taken to apply the accelerated method: 1) The required fidelity for the complete network Ω must be defined and Ω must be split into Ω low -and Ω high -regions, such that Ω high ∪ Ω low = Ω and Ω high is as large as possible.Here, Ω low is the set of regions that require a good fidelity method and should be simulated on a low abstraction level, whereas Ω high is the set of regions that can be simulated relatively accurately using methods on a high abstraction level.
2) The resulting pressure and velocity fields of the corresponding simulation methods must be equal (or at least in close vicinity) on the interface Γ = Ω high ∩ Ω low to ensure continuity of the complete solution φ in the complete network Ω.This means that the pressure and velocity values must be communicated between the simulation methods and subsequently updated in Ω low and Ω high .In the next section, the implementation details on both steps are illustrated.

IV. IMPLEMENTATION DETAILS
To properly describe the implementation details for the method proposed above, the continuous channel-based microfluidic device sketched in Fig. 1 is used as a running example.Furthermore, the LBM and MNA (as reviewed in Section II) are used as simulation methods for Ω low and Ω high , respectively.The idea proposed above can be realized as follows.

A. Step 1: Identifying the Required Fidelity
Recall that the first step aims at identifying (ideally many) Ω high -regions where a high abstraction simulation method is sufficient and, hence, can be utilized to accelerate the required computations.As sketched before in Fig. 1, those regions can usually be identified easily.For example, straight channels belong to Ω high , while, e.g., crossings and T-junctions better remain in Ω low .Γ should be located inside a straight channel, sufficiently far from, e.g., a crossing or junction, where the flow is organized and the pressure can be regarded as a point-value along the channel.The resulting separation of the network is illustrated in Figure 2a, where the two-dimensional regions represent Ω low , and the lines and nodes constitute Ω high .The network during the initial iteration; Ω low is replaced by fully connected graphs, and the resulting network is used to find q 0 .Please note that the identification step is not restricted to straight channels, crossings and T-junctions.Also, arbitrary channel shapes or components on the microfluidic device (e.g., heaters, mixers, droplet generators, etc.) could be identified as Ω low as long as corresponding high abstraction level simulation methods are available.

B. Step 2: Defining the Communication
The second step aims at ensuring that the complete solution φ is continuous in Ω.Hence, the correspondingly obtained values (here, pressure and velocity fields) from both simulation methods need to be adequately communicated from/to Ω high and Ω low , and they must be subsequently updated in the respective regions.This is similar to simulation methods of multiphysics problems [26], [27].Eventually, once these values align, we obtain a converged complete solution φ.In the proposed method, this is accomplished using an iterative method, i.e., a method that tries to find a local solution (fixed point) iteratively [28], [29].More precisely, where, in this case, q n is the quantity on Γ that is communicated between Ω low and Ω high at timestep n, and f (q) is the update function, given by the iterative method.The stability and convergence characteristics of the iterative method highly depend on the quality of the initial approximate solution, i.e., the initial condition q 0 .To find the initial condition q 0 , the complete network is first solved completely using the MNA, i.e., the high abstraction method.In Fig. 1, Ω low consists of crossings and T-junctions but could, in reality, contain any arbitrary shape that can be simulated using CFD.Therefore, the regions in Ω low are replaced by fully connected graphs, where each in-/outlet is treated as a node.This replacement is depicted in Fig. 2b.The resulting network can be solved using only the MNA.This is not the solution φ to the actual problem, since we abstracted Ω low by fully connected graphs.However, the solution to this abstracted problem can be used as an initial condition q 0 for the iterative method.
Having that, the communication can be performed in either of the two ways sketched in Fig. 3.A low abstraction solver (such as the LBM solver) calculates the pressure and velocity fields directly, whereas a high abstraction solver (such as the MNA solver) calculates the flow rate in a channel rather than the velocity field.This means that if we map the flow rate from the MNA solver to the LBM solver (Figure 3a), we need to extrapolate the flow profile based on the flow rate, whereas the reverse mapping (Fig. 3b) can be done directly (provided that Γ is located sufficiently far in a straight channel section, such that the pressure is sufficiently uniform on the channel cross-section and can be treated as a point value).Regions that belong to Ω high but are not connected to ground nodes (the green nodes in Fig. 2a) need to communicate according to Fig. 3a in at least one node.If this is not the case, the absolute pressure, i.e., the pressure difference with respect to the reference pressure p 0 at the ground nodes, is not propagated correctly.These regions are highlighted in red in Fig. 2a).To ensure a converged complete solution, we used an iterative method based on Successive Over-Relaxation (SOR, [28]) to find the values for pressure and velocity on Γ.From the initial condition q 0 , the boundary conditions of the LBM solver are updated in every iteration according to the newly found values from the MNA, corrected by a relaxation factor α, i.e., Here, q n high is evaluated using the MNA with the most recent pressure and velocity information obtained from the LBM solver.The LBM is in itself also an iterative solver, and the frequency at which the boundary conditions are updated influences the stability of the LBM.To ensure stability in the LBM, the values of the boundary conditions are only updated every θ "LBM iteration steps".The iterative approach in Eq. ( 8) is solved until the convergence criterion is met, where ϵ can be chosen arbitrarily small (until machine precision is reached).

V. CASE STUDIES AND EVALUATION RESULTS
The approach to accelerate CFD simulations, as proposed above, has been implemented for the considered continuous channel-based microfluidics.As corresponding simulation tools, we used OpenLB v1.5 [30] for the LBM and an in-house implementation of the MNA (which, as discussed above, are used as the low and high abstraction simulation methods, respectively).The source code of the proposed method, that has been developed for this work, is available online (available at: https://github.com/cda-tum/mmft-hybrid-simulator; accessed on 6 October 2023) .
Using the resulting implementation, several case studies were conducted to evaluate whether the idea proposed in this work indeed yields an improvement in CFD simulations.This section summarizes the respectively obtained findings.To this end, first the setup of the case studies is reviewed.Afterward, the obtained results are presented and discussed.

A. Setup: Considered Cases and Parameters
In our case studies, four different continuous channel-based microfluidic networks, as shown in Figure 4 (denoted Network 1-4 in the following), were considered.All networks are two-dimensional.For each of these networks, we considered different amounts of disconnected regions in Ω low , as well as different lengths of the channels that constitute Ω highproviding a proper variety of test cases with different coverages of Ω low -and Ω high -regions.More precisely, following the discussion from Section IV, all junctions and crossings of channels are identified as Ω low -region (indicated by the red dotted squares in Fig. 4), and the connecting channels are identified as Ω high .Γ was always placed inside a straight channel at a distance of two channel thicknesses from the corresponding junction or crossing.The length l of all channels is subsequently set to 1, 2, 3, and 4 mm for all four networks, whereas the Ω low -regions remain constant.Overall, this leads to a total of 16 separate cases.
The channels of all networks are rectangular with a width of 100 µm, and all networks have the inlets located on the left-hand side, with pressure boundary conditions of 1000 Pa, and outlets on the right-hand side, with pressure boundary conditions of 0 Pa.For all test cases, the fluid inside the network is an incompressible homogeneous fluid with a density of 1000 kg/m 3 and a kinematic viscosity of 1•10 −6 m 2 /s.
The test cases are solved on a regular grid with a resolution of 20 grid cells across the width of each channel, where applicable (i.e., the resolution of the CFD simulation and Ω low ).Additionally, for the proposed method, we set the values θ = 10, ϵ = 0.01, and α = 0.01 for Networks 1 and 2, and 0.003 for Networks 3 and 4. All simulations were performed without compiler optimization on a single CPU core (no parallelism) of an AMD Ryzen Threadripper PRO 5955WX CPU [31].

B. Obtained Results
In order to evaluate the performance of the proposed method, two aspects are essential: The runtime required for the respective CFD simulations (as we are aiming to accelerate them), as well as the accuracy (as potential accelerations ideally should yield the same results).
Concerning the former, Table I lists the respectively obtained results.Here, for all networks from Fig. 4 (listed from left to right), as well as for all channel lengths (listed in the rows), the correspondingly required runtimes of the original as well as the proposed method are provided.Additionally, the resulting speed-ups obtained by the proposed method are listed.As mentioned previously, the LBM was used in this work as a representative for CFD simulation methods.However, the computational complexity of the LBM is similar to that of, e.g., the FVM [20] and speed-ups of similar order of magnitude can be expected for other simulation methods.
Concerning accuracy, direct comparisons of the pressure and velocity fields obtained with CFD simulations and the proposed method are given in Figures 5-8.Fig. 5 shows all the obtained results for Region 1a of Network 1 (Fig. 4a) for all four channel lengths l, i.e., the complete Network 1-column in Table I.Figs. 6 to 8, respectively, show the obtained results for all Ω low -regions in Networks 2, 3, and 4 Fig. 4b-d, at channel length l = 1, i.e., the top row in Table I.Finally, we list the pressure values and velocity magnitudes obtained by both approaches for all test cases in Tables II and III.These values were taken in the Ω low -regions (labeled 1a, 2a, ..., 4c, 4d in Fig. 4) at the measuring points indicated by the black dots in Fig. 4. Since the networks were simulated for four different channel lengths l, each measuring point has four pressure values.The results clearly confirm the improvement and benefit of the proposed acceleration method.First, we can see from the numbers summarized in Tables II and III that both the original CFD simulation, as well as the proposed method, more or less provide the same results.That is, using higher levels of abstractions for Ω high -regions does not significantly affect the simulation results.At the same time, this enables impressive speed-ups.In fact, the numbers summarized in Table I do not only show that the proposed method always generates the results faster than the original CFD method, but also that it constantly does so by several factors-in many cases even by up to three orders of magnitudes.

C. Discussion
The obtained results, as summarized above, clearly show the benefits of the proposed method.On top of that, they also provide further, more detailed insights, as well as implications and ideas for further extensions.These are discussed in the following.
Firstly, with respect to the accuracy, the results presented in Figs. 5 to 8 and Tables II and III show a strong alignment between the proposed method and the corresponding CFD simulations.From Fig. 5, it can be concluded that there seems to be a slight improvement in the pressure field as the channel length l increases, but this is negligible.The results in Figs. 6 to 8 show that the method is applicable for problems with multiple disconnected Ω high -regions, with a negligible increase in inaccuracy for Ω high -regions that are located between others, such as 3b, 4b and 4c.Following from the pressure contour lines obtained from the CFD approach, which generally appear straight towards Γ, the location of Γ can be said to be sufficiently far away from the crossings and junctions and did not strongly influence the performance of the proposed method.The pressure contour lines obtained from the proposed approach, however, consistently show a curvature near Γ, which can be attributed to numerical intricacies of the underlying implementation of the boundary conditions.This slight inaccuracy can also explain the propagated error to Regions 3b, 4b, and 4c.
The exploitation of higher levels of abstraction is possible due to the availability of reduced order modeling methods for microfluidic flow.In this work, the Hagen-Poiseuille law was used to represent the fluid flow through straight two-dimensional channels.Using a different method of high abstraction, that models flow through three-dimensional channels, such as presented in [32], would allow the proposed method to be extended to three dimensional problems.With this added dimension and, therefore, even worse computational complexity of CFD simulations, the speed-up of the proposed method can be expected to be of even higher orders of magnitude.Provided that methods of high abstraction exist, the proposed method can even be extended to include other physical phenomena, such as diffusion, heat dissipation, or droplets.For straight channels, diffusion and heat dissipation can simply be modeled as a time-dependent transport across the channel width, and high abstraction simulation approaches for droplets are performed by adding the droplet hydraulic resistance to the channel [33].Further work is required to include additional physical phenomena in the proposed method.

VI. CONCLUSION
In this work, we proposed an accelerated CFD simulation method for microfluidic devices.The core idea was to utilize higher levels of abstraction whenever possible to improve the simulation runtime while maintaining the fidelity.We developed and implemented a prototype of the resulting simulation approach, using continuous channel-based microfluidic devices as a representative platform.Results obtained from corresponding case studies confirmed the promises of the proposed approach: using higher levels of abstractions for the simulation did not significantly affect the fidelity of the simulation results, but allowed for substantial speed-ups of up to three orders of magnitude.Based on this premise, and with the inclusion of other physical phenomena such as diffusion, heat dissipation, or multiphase flow (e.g., droplets), similar accelerations can be expected for further microfluidic platforms, which is left for future work.Overall, this provides the foundation for more research toward exploiting higher levels of abstraction for simulating microfluidic devices in future work.

Fig. 1 :
Fig. 1: Example network for continuous channel-based microfluidics.The pressure p and velocity fields u are shown in detail for a crossing (red) and a straight channel section (green).

Fig. 2 :
Fig. 2: The example network in the proposed method.(a) The separation of the example network into Ω low and Ω high .(b)The network during the initial iteration; Ω low is replaced by fully connected graphs, and the resulting network is used to find q 0 .

Fig. 3 :
Fig. 3: Communication schemes for the pressure and flow fields.(a) Communicate the flow rate from Ω high to Ω low and the pressure vice versa.Here, the flow field information must be extrapolated from the communicated flow rate.(b) Communicate the flow field from Ω low to Ω high and the pressure vice versa.

Fig. 4 :
Fig. 4: The networks of the considered case studies.The number of separate regions in Ω low increments with the case studies, starting at one region (a) for Network 1 and ending at four regions (a-d) for Network 4.

TABLE I :
Required runtimes of the original CFD simulations and the proposed method with corresponding speed-ups.

TABLE II :
Obtained pressure values of all test-cases at the measuring points (as denoted in Figs.4a to 4d) obtained from the original CFD simulation and the proposed method.