Inverse design of microchannel ﬂuid ﬂow networks using Turing pattern dehomogenization

Microchannel reactors are critical in biological plus energy-related applications and require meticulous design of hundreds-to-thousands of fluid flow channels. Such systems commonly comprise intricate space-filling microstructures to control the fluid flow distribution for the reaction process. Traditional flow channel design schemes are intuition-based or utilize analytical rule-based optimization strategies that are oversimplified for large-scale domains of arbitrary geometry. Here, a gradient-based optimization method is proposed, where effective porous media and fluid velocity vector design information is exploited and linked to explicit microchannel parameterizations. Reaction-diffusion equations are then utilized to generate space-filling Turing pattern microchannel flow structures from the porous media field. With this computationally efficient and broadly applicable technique, precise control of fluid flow distribution is demonstrated across large numbers (on the order of hundreds) of microchannels.


Introduction
Microchannel flow structures are found in a range of important industrial applications involving water purification (Wang et al. 2014), pharmaceuticals (Gutmann et al. 2015), electronics (Chen and Cheng 2002), and green chemistry (Lerou et al. 2010). Such systems require careful handling of fluid species to control reaction processes. Particularly, the flow distribution of fluids is critical, and intricate space-filling channel structures are employed for flow control through single-or multi-layered (e.g., 100-10,000) microchannels. Design optimization schemes for uniform fluid delivery to arrays of microchannels include approximation techniques, heuristic strategies, bifurcation methods, and gradient-based algorithms; see Rebrov et al. (2011). One approximation technique (Commenge et al. 2002) treats flow friction via a resistive network of ducts to understand pressure drop across a flow distribution chamber (or manifold). Using analytical expressions, the influence of the manifold geometry is understood, and flow uniformity is optimized for simplified (e.g., trapezoidal) geometries. A heuristic technique (Luo et al. 2015) involves perforated baffles upstream of microchannels to homogenize fluid distribution, although added pressure drop is a concern. The use of constructal theory (Bejan and Errera 1997) is popular with some (Senn and Poulikakos 2004) using tree-like flow distributors for fuel cells and others (Tondeur and Luo 2004), applying fractal manifolds to hundreds of channels in an adsorbent monolith. An associated challenge may be the spatial requirements for series configured branching. Borrvall and Petersson (2003) pioneered gradient-based topology optimization for Stokes flow, and this was extended to higher Reynolds numbers (Gersborg-Hansen et al. 2005;Olesen et al. 2006). The method was employed (Dede et al. 2014) to optimize manifolds with less than 10 fluid outlets in an electronics heat sink. Recently, researchers (Liu et al. 2011;Zeng and Lee 2018) incorporated mass flow rate constraints for a slightly greater number of outlets, e.g., ∼ 20. We leverage these formulations, where Navier-Stokes laminar incompressible fluid flow (for Reynolds numbers, Re < 2100) in an idealized porous medium is assumed with a friction force proportional to the fluid velocity, viz. Darcy's law, (Olesen et al. 2006;Kaviany 1995), Fluid incompressibility is governed by (1), and flow in the idealized porous medium is defined by (2). The fluid pressure and velocity vector state variables are given by P and u, respectively. The fluid dynamic viscosity is η. The effective inverse permeability is a function of a design variable, γ , and in standard formulations, Olesen et al. (2006), is interpolated using the convex function α = α s , A low permeability quasi-solid state exists for γ → 0 by α max = 1/l 2 Da , where Da is the Darcy number, and l is the characteristic length. A fluid state, where the friction force term goes to zero when γ → 1, emerges with α min = 0. The goal in adopting (3) with the tuning parameter, q, is to obtain completely fluid or solid states. However, challenges arise in selecting the appropriate Darcy number, Da, to eliminate gray-scale designs that permit undesirable flow seepage through quasi-solid material.
For a design space, , an accepted (Borrvall and Petersson 2003;Olesen et al. 2006;Gersborg-Hansen et al. 2005;Dede et al. 2014) objective function, f o , is to minimize power dissipation or flow resistance, as follows: A volume constraint is typical to control the amount of fluid in the result (Olesen et al. 2006;Gersborg-Hansen et al. 2005;Dede et al. 2014). However, spacefilling structures are not obtainable since channel-to-wall spacing is not enforced. In contrast, removal of the volume constraint produces large open flow channel designs to naturally minimize flow resistance.
For uniform mass flow through discrete fluid outlets, a mass flow rate constraint is common (Liu et al. 2011;Zeng and Lee 2018): whereṁ k andṁ t,k are the individual and target mass flow rate, respectively, at the kth outlet, and δ is an outlet flow rate error tolerance. To represent large gradients in the flow solution (e.g., near outlets), many channels require an extremely fine computational grid, per Reddy and Gartling (2000). Thus, effective optimization strategies for fluid flow control to hundreds of microchannels in application is still a critical field of research. Complex space-filling microchannel patterns are common, and for arbitrary domains, configurations may not be expeditiously found using approximation, rule-based, or modern gradientbased techniques. Complex patterns also exist in nature at multiple scales (e.g., mammalian markings, fish skin, and seashells), and reaction-diffusion equations replicate irregular spatiotemporal Turing patterns (Pearson 1993;Gray and Scott 1985;Garikipati 2017;Kondo and Miura 2010). Here, we combine gradient-based porous media optimization in an arbitrary fluid flow domain with a reaction-diffusion model for computationally efficient development of intricate space-filling microchannel architectures. Explicit modeling of channels is abandoned, and channel synthesis is realized through a Turing pattern generation algorithm.

Methodology
The approach comprises two steps. First, to scale up to manifolds with hundreds of outlets, the problem is re-framed to design a homogenized porous fluid flow structure, where all material states are physically feasible. Following Kim and Kim (1999), the porous media is parameterized in two-dimensions (2-D) based on a spatially varying local microchannel structure; see Fig. 1. The porosity, , and permeability, κ, are as follows: where w c and w w are the channel and wall widths, respectively. A linear interpolation function for the channel width relates any porous media state to a microstructure, as follows: where w c min and w c max are minimum and maximum microchannel widths, respectively. Combining (6) and (7), and assuming w w is constant, a new inverse permeability expression, α → α n , is as follows: Thus, 0-1 (fluid solid) designs are not required, flow seepage is not a concern, a volume constraint is unnecessary, and hundreds of discrete fluid flow outlets are exchanged for an aggregated fluid outlet with specified flow profile. The second step of the algorithm to find the Turing pattern is an expansion of the anisotropic thermal-composite design approach explained in Petrovic et al. (2018). The Turing reaction-diffusion system is a mathematical model of the morphogenesis of the embryo proposed by Turing (1952), involving two interacting hypothetical chemical substances U and V, which diffuse in the space around and enhance or suppress the reproduction of themselves (Kondo and Miura 2010). The dimensionless equations for this process are as follows: where R u (U, V ) and R v (U, V ) are interactive reaction terms. In this study, the reaction terms are augmented following (Kondo and Miura 2010;Petrovic et al. 2018): For 2-D flow systems, we extend the diffusion coefficients as 2 × 2 tensors, D u and D v , with anisotropic diffusion terms derived from local permeability field values and perturbed over time between weakly and strongly anisotropic states. By aligning the principal axis of the diffusion tensors with the fluid velocity vector, the microstructure pitch (i.e., the distance between microstructures) and length are controlled with length periodically elongated along the fluid flow direction.
For example, D u is written using the normalized fluid flow velocity vector,ū, as follows: where ⊗ is the dyadic product operator. The coefficients in (15) are defined as L u = l 2 u W u and W u = (w u w) 2 with l u set to control the magnitude of anisotropy and w u specified as a constant value that associates the channel pitch and generated microstructure pattern. By specifying the channel pitch, w, the lateral magnitude of the diffusion, i.e., W u , is proportional to w 2 , and thus, we recover the porous media permeability distribution, w = w c + w w in (6), where w w is constant in the numerical examples. Similarly, D v is defined as follows: where L v = l 2 v W v and W v = (w v w) 2 . The magnitude of the Turing pattern anisotropy repetitively varies temporally by changing l u between ∼ 1 and 10.
Larger anisotropy generates a pattern strongly affected by the fluid velocity vector field, and smaller anisotropy generates a new more tightly packed pattern based on the previous final state; refer to the Supplementary Movie. Figures 2 ae show an example with an optimized fluid flow velocity field in an effective porous medium with varying permeability and the evolution of the Turing pattern microstructure. The algorithm beneficially combines homogenization of the design space for the computationally expensive flow field optimization with efficient dehomogenization of the porous media. This dehomogenization process may uniquely vary, and other approaches to dehomogenization for topology optimization in linear elasticity are available; see for example (Groen and Sigmund 2018;Allaire et al. 2019). However, here, to render the structure on to a physical scale, a self-organizing system based on Turing patterns is selected.

Implementation
The algorithm is implemented in COMSOL v.5.3a. The optimization objective function is set to minimize an equally weighted linear combination of two terms including (4) and variation of the fluid velocity normal to the outlet boundary.
The latter term is derived from the uniform mass flow constraint, (5), as follows: whereū j is the average velocity in the j -direction, and o is the outlet boundary. Channel and wall width values, more precisely, the minimum and maximum value of the channel width and the wall width, are specified to parameterize the effective inverse permeability, α n , of the flow space. The flow optimization uses a method of moving asymptotes (MMA) optimizer that converges within ∼ 100 iterations. From the porous media field results, the anisotropic diffusion coefficient tensors for the reactiondiffusion equations are determined, and the equations are propagated through time (∼ 1800 s) to generate the Turing pattern microstructure.

Results
A 200 × 100 mm 2 2-D space is considered in the upper image of Fig. 3a in a first numerical experiment with air flow at 20 • C. A uniform +y-direction fluid velocity of 0.2 m/s (Re = 133) is fixed over the 10-mm wide fluid inlet  positioned 15 mm from the domain lower left corner. This asymmetric inlet-to-outlet configuration represents typical microchannel reactors and produces significant outlet flow maldistribution; see Commenge et al. (2002). A zero pressure outlet boundary condition (BC) is applied along the top edge of the flow space. A Dirichlet BC for the reaction-diffusion model is applied along the same top edge of the domain and on a horizontal inner boundary slightly (3 mm) below the top edge (and meshed to have element edges at the desired segment) to enforce a precise channel width distribution. A second Dirichlet BC is applied on the remaining domain boundaries. The minimum and maximum channel width is w c min = 0.6 mm and w c max = 1.8 mm, respectively. The wall width between channels is uniformly fixed to w w = 0.6 mm. Thus, porosity, , ranges from 0.5 to 0.75 throughout the domain, while permeability, κ, ranges from 1.5E-8 m 2 to 2.025E-7 m 2 . The channel width at the top outlet boundary is set to 0.6 mm, and ∼ 166 uniform width outlets are implicitly defined. The parameter values for the Turing pattern generation are provided in Table 1.
The gray-scale fluid flow optimization result for the effective porous media is shown on top in Fig. 3a including streamlines and normalized pressure contours with the Turing pattern microstructure shown below. Note that intricate patterned island wall structures are built following the above parameters. Fluid velocity contours from a flow verification analysis under the same inlet/outlet BCs are shown in Fig. 3b. Observe that the average variation, in the mass flow at each outlet,ṁ k , is very low and within 1.3% of the mass flow average,ṁ avg , across the n = 166 outlets; the maximum variation, max[| (ṁ k − m avg )/ṁ avg |], at a single outlet is 4.6%. For comparison, the maximum variation inṁ k for the flow domain without any microstructures is 52.9%; see Fig. 3b. Additionally, the mass flow performance of the homogenized grayscale design across the 166 outlets is shown in Fig. 3b Table 1 Turing pattern generation parameters with very good agreement between the homogenized and dehomogenized designs.
To compare the computational effort involved in our approach with conventional flow channel topology optimization, cf. Olesen et al. (2006), Liu et al. (2011), and Zeng and Lee (2018), a simple 100 × 100 mm 2 2-D flow space is considered in a second numerical experiment. Here, the 10-mm wide fluid inlet is positioned at the center of the lower edge of the domain. A uniform +y-direction fluid inlet velocity of 0.2 m/s (Re = 133) is again assumed; the remaining BCs follow the prior example. The same microchannel and wall widths are further assumed for 83 implicitly defined outlets; the optimized Turing pattern microstructure is shown in Fig. 4. Three flow structures from conventional topology optimization (TO1, TO2, and TO3) are also shown with 10, 20, and 30 explicitly defined fluid outlets obtained using (5). Convergence of designs requires 70 iterations using Da =5E-7 in (3). To relax the optimization problem and avoid local solutions, q is set to 0.01 and ramped up to 0.05 after the 30th iteration in 2E-3 increments. Also, design variable regularization follows (Kawamoto et al. 2011) with a Heaviside bandwidth, H d , initially set to 0.6 and ramped down to 0.2 after the 20th iteration in 1.33E-2 steps. For explicitly defined outlets, the computational cost measured in CPU hours using a laptop with 12, 2.71 GHz processors grows quadratically in Fig. 4 due to the progressively finer mesh requirement at the outlets. However, the total computational cost required to optimize the porous media field plus generate the microstructure Turing pattern for a design of any number of implicitly defined outlets remains constant, assuming an initially sufficiently refined mesh to ensure solution accuracy in both algorithm solution segments. Thus, our computational procedure is efficient and additionally bounds the minimum and maximum channel widths throughout the flow domain. Table 2 provides details on the average outlet mass flow variation, maximum outlet mass flow variation, and pressure drop from an additional flow verification analysis for each design. The computational cost to obtain the various designs is also noted. Note that all pressure drop, P , values are normalized by the lowest value for the 30 outlet design. For the fully post-processed designs optimized using conventional topology optimization, the average outlet mass flow variation ranges from 7 to 12%, while the maximum outlet mass flow variation is ∼ 20%. For the post-processed designs obtained using the dehomogenization procedure, the average outlet mass flow variation is 0.75-2%, and the maximum outlet mass flow variation is less than 5.9%. Observe that the pressure drop for the dehomogenized Turing pattern designs is approximately 4-6× that of the conventional topology optimized designs. This greater flow resistance is logically related to the assumption of porous media throughout the Fig. 4 Computational cost versus number of outlets for the full algorithm solution for a 100 × 100 mm 2 manifold. TO, conventional topology optimization algorithm; RD, reaction-diffusion algorithm. Light regions = fluid, black regions = solid design space combined with an eventually large number of very small diameter outlet channels as opposed to fewer flow channels with unconstrained maximum size. Here, w c max and w c min are fixed in all cases, and the 83 outlet case shown in Fig. 4 is representative of all outlet configurations with the exception of the microstructure as it approaches the outlets. However, as a future extension of this work, we may also modify/optimize w c max , w c min based on further space-filling microstructure, fabrication, and flow resistance constraints.

Conclusions
An optimization method was proposed to design spacefilling Turing patterned microstructures for fluid flow control through a porous media field. Reaction-diffusion equations form the basis of the microchannel dehomogenization post-processing technique, where the porous media field is linked to explicit representation of the rendered microstructure. Numerical experiments highlight the capability to produce precise flow control for manifolds involving hundreds of microchannel outlets with a trade-off to higher pressure drop due to the space-filling microstructures. The method is effective, flexible, and may be applied to large arbitrary geometries. Other porous media parameterizations are feasible, including permeability mappings that account for the flow structure out-of-plane depth, opening opportunities for design in three dimensions. Relevant uses of this methodology include the design and additive fabrication of microchannel reactors, which are prevalent across biological and energy-related applications.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.

Replication of results
Refer to the Supplementary Movie as a detailed intermediate result. Due to institutional constraints, the source code is unavailable. However, further algorithm details are available upon request from the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.