Isogeometric Hierarchical Model Reduction for advection-diffusion process simulation in microchannels

In the last decade, microfluidics has gained increased interest, becoming a key technology in the fields of biology and medical research, with various applications such as biosensing, diagnostics, drug discovery, regenerative medicine, tissue engineering (Nguyen et al. [2019], Sackmann et al. [2014], Bhatia and Ingber [2014], Davis [2008]). The strength of microfluidics relies on the possibility to substantially reduce the sample volume by using disposable miniaturized devices which allow to replace large-scale conventional laboratory instrumentation, thus reducing hardware costs, assuring low reagent consumption and high-speed analysis, while allowing huge parallelism.

In a microfluidic system, mixing of different transported passive species is a crucial issue for the optimization of chemical and biochemical reactions involved in the process of medical diagnostics, drug discovery, chemistry production and proteomics (Suh and Kang [2010]). However, at the microscale, mixing is often difficult to be achieved in common practice, since microfluidics is characterized mainly by very low Reynolds numbers (ranging from less than unity up to a few hundreds), and cannot take advantage of turbulence to improve mixing efficiency. Under laminar flow, mixing is mainly a diffusion-driven phenomenon characterized by long time scales and high diffusive lengths, which can become prohibitive in view of a numerical modeling.

To overcome these limitations, great efforts have been made to develop new techniques in order to achieve rapid and efficient laminar flow mixing in microsystems (Hardt et al. [2005]). In particular, chaotic advection can be exploited to enhance the diffusive process in microfluidics. Chaotic advection is obtained by stretching and folding the interface between miscible streams, thus reducing the diffusive length, and improving the mixing efficiency, also at very low Reynolds numbers. In passive micromixers, chaotic advection is usually achieved by properly designing the three-dimensional geometry of microchannels (Liu et al. [2000], Lin et al. [2007]). It is well known that in a curved channel the flow undergoes a centrifugal displacement of the maximal axial velocity (Dean [1928]), leading to the onset of secondary flows, which in turn enhance mixing. Thus, several micromixer designs are based on curved microchannels, such as spiral microchannels (Sudarsan and Ugaz [2006]), U-shaped microchannels (Gigras and Pushpavanam [2008]), or clothoid-based geometries (Pennella et al. [2012]).

Following this principle, in this study an S-shaped microchannel with fixed curvature is considered, in the simplest case of a single passive species. In such a way, we recover the conventional setting to apply a Hierarchical Model (HiMod) reduction. As a matter of fact, HiMod proved to be an effective tool (Guzzetti et al. [2018], Brandes Costa Barbosa and Perotto [2020]) to model phenomena exhibiting a main dynamics (here represented by the advection of the passive scalar along the microchannel), in the presence of local secondary dynamics evolving along the transverse sections (in the case of interest, represented by the mixing effects induced by the geometry). Analogously to other model reduction procedures (we refer, e.g., to Chinesta et al. [2014], González et al. [2010] and to Perotto et al. [2020] for a comparison between HiMod and another well-established model reduction technique), a HiMod reduction exploits a standard separation of variables and describes the mainstream and the secondary dynamics with different approximation schemes. In the first proposal, the mainstream is discretized by one-dimensional finite elements, while the transverse dynamics are modeled by a modal expansion (Ern et al. [2008], Perotto et al. [2010], Perotto and Zilio [2013], Perotto [2014]). For the specific context at hand, we resort to the variant of the original approach which employs an isogeometric discretization along the leading direction, in order to deal with generic geometries (Perotto et al. [2017], Brandes Costa Barbosa and Perotto [2020]).

Independently of the selected discretizations, HiMod reduction considerably contains the computational effort, in particular when compared with full three-dimensional models, without quitting accuracy. Indeed, a HiMod expansion leads to commute the full model into a system of coupled one-dimensional problems solved along the leading direction, whose coefficients include the effect of the transverse dynamics. This ensures several computational simplifications, especially in the presence of complex geometries, primarily the computational domain we have to discretize which is a one-dimensional instead of a three-dimensional one. In terms of accuracy, the HiMod approximation can be arbitrarily enriched by properly increasing the number of modal functions used to describe the secondary transverse dynamics (Perotto et al. [2010], Perotto and Veneziani [2014], Perotto and Zilio [2015], Aletti et al. [2018]).

The paper is organized as follows. Section 2 introduces the partial differential equation model used to describe the transport of the tracked species inside the S-shaped microchannel. Section 3 is devoted to set up the isogeometric HiMod reduction, by addressing (i) the geometric characterization of the computational domain, which distinguishes between the mainstream and the secondary directions; (ii) the definition of thefp discrete HiMod space; (iii) the HiMod algebraic formulation. Numerical results are discussed in Section 4, both in terms of accuracy and computational saving. Finally, Section 5 summarizes conclusions and possible future perspectives.

The problem we consider models the transport of a passive scalar in the S-shaped microchannel with two fixed curvatures, $\Omega\subset\mathbb{R}^{3}$, sketched in Figure 1. To this aim, we resort to a standard advection-diffusion (AD) equation, that is formulated as⁵⁵5Throughout the paper, standard notation are adopted for the function spaces (Ern and Guermond [2004]).

where $u$ denotes the concentration of the passive scalar, the convective field $\mathbf{b}\in[L^{\infty}(\Omega)]^{3}$, with $\nabla\cdot\mathbf{b}\in L^{2}(\Omega)$, corresponds to the velocity of the fluid in the microchannel, $f\in L^{2}(\Omega)$ models the action of an external force, $\nu\in\mathbb{R}^{+}$ is the diffusion coefficient of the transported species, $u_{in}\in H^{1/2}(\Gamma_{in})$ is the concentration at the inlet $\Gamma_{in}$, with $\mathbf{z}$ the vector collecting the spatial coordinates in $\Omega$, $\mathbf{n}$ the outer unit normal to the boundary $\partial\Omega$, $\Gamma_{in}$, $\Gamma_{out}$ the inlet and the outlet sections, $\Gamma_{w}$ the lateral boundary, so that $\partial\Omega=\Gamma_{in}\cup\Gamma_{w}\cup\Gamma_{out}$.
We observe that, since the dynamics in a microchannel is usually characterized by a very small Reynolds number, we can assume that the fluid is incompressible, namely $\nabla\cdot\mathbf{b}(\mathbf{z})=0$ in $\Omega$, and that the fluid flow is laminar, so that, in practice, we can compute $\mathbf{b}$ by solving a steady Navier-Stokes problem (see Section 4).

In view of the application of the Hierarchical Model (HiMod) reduction technique, we introduce the weak formulation of problem (1). To this end, we introduced the Sobolev spaces $H^{1}(\Omega)=\{v\in L^{2}(\Omega):\nabla v\in\ [L^{2}(\Omega)]^{3}\}$ and $H^{1}_{0,\Gamma_{D}}(\Omega)=\{v\in H^{1}(\Omega):v|_{\Gamma_{D}}=0\}$, with $\Gamma_{D}=\Gamma_{in}\cup\Gamma_{w}$. Thus, the weak formulation of the reference AD problem reads: Find $u\in V$, with $V=H^{1}_{0,\Gamma_{D}}(\Omega)$, such that

where $a:V\times V\longrightarrow\mathbb{R}$ and $F:V\longrightarrow\mathbb{R}$ are the bilinear and the linear forms associated with model (1), namely

with $R_{u_{in}}$ a suitable lifting of the boundary data $u_{in}$ on $\Gamma_{in}$.
Thanks to the assumptions introduced on the problem data, we can guarantee the well-posedness of the weak problem (2) by applying the Lax-Milgram lemma (Ern and Guermond [2004]).

Hereafter, we will refer to (2) as to the full problem, the HiMod reduction will be applied to.

In this section, we apply the HiMod reduction to the full problem (2), by discussing the three main steps typical of a HiMod formulation: (i) the geometric characterization of the computational domain that distinguishes between a leading and a transverse direction (Section 3.1); (ii) the definition of the function spaces associated with either directions and their combination in the setting of the HiMod discrete space (Section 3.2); (iii) the resulting algebraic formulation, which reduces to a system of one-dimensional (1D) coupled problems (Section 3.3). For further details, we refer the reader, e.g., to Ern et al. [2008], Perotto et al. [2010], Perotto [2014].

In this section we apply the HiMod discretization to model the transport of a passive scalar in an S-shaped microchannel, by hierarchically reducing problem (1) in the computational domain $\Omega$ displayed in Figure 1.

The geometry of the microchannel is characterized by a length L of $1410$ $\mu$m, with a circular section of diameter $D=100$ $\mu$m, a first radius of curvature $R_{1}$ equal to $455$ $\mu$m and a second radius of curvature $R_{2}$ equal to $204$ $\mu$m.

Concerning the problem data, we assume that the external force $f$ in (1) is null, while the convective field $\mathbf{b}$ is obtained by solving the steady-state Navier-Stokes equations in $\Omega$. To this aim, we use the general purpose CFD code Fluent (ANSYS Inc., Canonsburg, PA, www.ansys.com), based on a finite volume discretization. The Navier-Stokes simulation is completed with the following boundary data: we prescribe a constant velocity of $2$ $m/s$ at the inlet section $\Gamma_{in}$; a traction-free condition at the outlet section $\Gamma_{out}$; a no-slip condition on $\Gamma_{w}$, thus assuming the microchannel wall to be rigid. The dynamic viscosity of the fluid is set to $\mu=10^{-3}$ $kg/ms$, while the density is $\rho=998$ $kg/m^{3}$. Therefore, the resulting regime for the fluid flow is laminar, as the Reynolds number at the inlet section is approximately equal to $200$.

The diffusion coefficient of the transported species in (1) is $\nu=10^{-6}$ $m^{2}/s$, so that the Péclet number at the inlet section corresponds to ${\mathbb{P}}{\mbox{e}}=200$. Finally, as for the boundary conditions of the AD equation, the scalar distribution $u_{in}$ corresponds to a parabolic profile on the inlet section $\Gamma_{in}$, while we assign a homogeneous Dirichlet boundary condition on the microchannel wall $\Gamma_{w}$, and a homogeneous Neumann data on the outlet section $\Gamma_{out}$.

Concerning the geometric characterization of the computational domain in view of a HiMod reduction, the leading direction is aligned to the centerline of the microchannel (i.e., to the direction the convective field flows along), while the transverse sections coincide with the uniform circular cross-sections of the channel.
The NURBS description of volume $\Omega$ is built as follows. For each point $\mathbf{P}=(x_{p},\textbf{y}_{p})$ along the centerline (i.e., such that $x_{p}\in\Omega_{1D}$), we provide the function $\mathbf{p}:[0,2\pi]\to\mathbb{R}^{3}$ such that the curve

identifies the boundary of the transverse fiber, $\gamma_{x_{P}}$, of $\Omega$ centered at $\mathbf{P}$, where vectors $\mathbf{B}$ and $\mathbf{N}$ correspond to the binormal and to the normal unit vectors of the Frenet frame at $\mathbf{P}$, respectively. Such vectors can be easily computed by using, e.g., vmtk (www.vmtk.org).
Now, we focus on maps $\Phi$ and $\Psi$, relating the physical with the reference domain. First, we introduce a NURBS map between each side of the reference transverse fiber $\hat{\gamma}=(0,1)^{2}$ and a portion of the wall of the generic circular cross section. To this aim, the curve defined by (20) is built by varying angle $\theta$ in the four disjoint intervals $[-\pi/4,\pi/4)$, $[\pi/4,3/4\pi)$, $[3/4\pi,5/4\pi)$, $[5/4\pi,7/4\pi)$. Successively, each of the four resulting curves is mapped to one side of the unit square $\hat{\gamma}$, according to the correspondence highlighted in Figure 2. Afterwards, a NURBS parameterization of the internal part of the transverse fiber is carried out by means of a bilinearly blended Coons patch (i.e., via a suitable bilinear interpolation of the four arcs in Figure 2, see, e.g., Farin and Hansford [1999]), as implemented in the Matlab NURBS toolbox. Finally, maps $\Phi$ and $\Psi$ are obtained by repeating the procedure above for each point $\mathbf{P}$ along the centerline, as sketched in Figure 3.

To identify the discrete HiMod space in (7), we resort to a modal basis consisting of $64$ functions which describe the dynamics parallel to the transverse direction, while we discretize the flow aligned with the centerline of the microchannel by means of quadratic NURBS basis functions, characterized by a $C^{1}$ inter-element smoothness, after introducing $202$ uniformly-spaced knots along the reference supporting fiber $\hat{\Omega}_{1D}$. Panel (b) in Figure 4 shows the HiMod solution. Four cross sections of interest are highlighted: the inlet (#1), a section halfway through the first straight part of the microchannel (#2), a section at the beginning of the region with curvature $R_{1}$ (#3), and finally a section halfway through the region with curvature $R_{1}$ (#4). Furthermore, a cut in the plane which contains the centerline is also shown.

For comparison purposes, we run a simulation of the same AD problem with Fluent on a uniform unstructured computational grid of $\Omega$ consisting of $1.384.925$ tetrahedral elements and $273.421$ vertices. In particular, Fluent solves the Navier-Stokes and the AD equations concurrently. We adopt a second-order discretization scheme for the pressure, and a second-order upwind method to approximate both the velocity and the passive scalar concentration. The resulting solution is considered as the high fidelity reference solution, and is employed to establish the accuracy of the HiMod approximation.

A qualitative comparison between the HiMod and the reference solution is displayed in Figure 4. We observe a good qualitative agreement between the two solutions. We highlight that the data assigned at the inlet of the HiMod reduced model does not exactly coincide with the parabolic profile $u_{in}$ used in the reference context. Indeed, as discussed in Remark 2, the lifting $R_{u_{in}}$ in (4) is built upon the expansion of the Dirichlet data $u_{in}$ in terms of the modal and of the NURBS bases. This justifies the minor discrepancy in the concentration on the inlet section (#1). Such a slight mismatch is unavoidably propagated along the channel and explains the small error present on the downstream sections. In particular, the HiMod solution slightly overestimates the high fidelity concentration on section #2, while underestimation occurs on sections #3 and #4.

A more quantitative comparison is performed in Figure 5, which displays the pointwise distribution of the relative modeling error between the high fidelity and the reduced solution with respect to the $L^{2}(\Omega)$-norm, in correspondence with the plane and the four sections considered in Figure 4. We can infer that a small relative error (less than $1\%$) is achieved, which is often acceptable for the type of engineering applications analyzed in this chapter.
Finally, from a computational viewpoint, the gain provided by the HiMod approximation is confirmed by the number of the degrees of freedom (dofs) characterizing the two solutions, namely $1.384.925$ dofs for the high fidelity model to be compared with $12.928$ dofs for the reduced model.

An S-shaped microchannel used to enhance mixing by yielding chaotic advection in passive micromixers provides the ideal environment where to verify the computational performances of a HiMod reduction procedure, both in terms of reliability and computational efficiency. Indeed, it is common to distinguish in the microchannel geometry a leading direction, aligned with the channel centerline, and an orthogonal transverse direction, parallel to the cross section, coherently with the separation of variables exploited in the definition of the HiMod reduced space.
The intrinsic capability of Himod reduction to decouple the dominant dynamic from the secondary one through a differentiated discretization of the two directions leads to the resolution of a system characterized by a considerably lower number of unknowns, without quitting the accuracy of the reduced solution. In particular, for a $99.07\%$ reduction in the number of dofs with respect to a high fidelity finite volume simulation, a very good accuracy is exhibited by the HiMod solution, with a less than $1\%$ relative modeling error.

In a future perspective, the proposed approach may be extended to more complex problems in the microfluidics field, such as the study of AD problems under unsteady conditions or in the presence of more complex 3D geometries. Moreover, a model reduction approach could speed up the microchannel design phase, especially when different scenarios have to be investigated in order to optimize the mixing efficiency of the microsystem. The parametric version of the HiMod reduction could be instrumental for such a goal (Baroli et al. [2017], Zancanaro et al. [2021], Lupo Pasini and Perotto [2022]).

All the authors acknowledge the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Actions, grant agreement 872442 (ARIA, Accurate Roms for Industrial Applications). Moreover, SP thanks the PRIN research grant n.20204LN5N5 (Advanced Polyhedral Discretisations of Heterogeneous PDEs for Multiphysics Problems); UM thanks the MIUR FISR FISR2019_ 03221 CECOMES project; FB thanks the project “Reduced order modelling for numerical simulation of partial differential equations” funded by Università Cattolica del Sacro Cuore.