Design Applicable 3D Microfluidic Functional Units Using 2D Topology Optimization with Length Scale Constraints

Due to the limits of computational time and computer memory, topology optimization problems involving fluidic flow frequently use simplified 2D models. Extruded versions of the 2D optimized results typically comprise the 3D designs to be fabricated. In practice, the depth of the fabricated flow channels is finite; the limited flow depth together with the no-slip condition potentially make the fluidic performance of the 3D model very different from that of the simplified 2D model. This discrepancy significantly limits the usefulness of performing topology optimization involving fluidic flow in 2D—at least if special care is not taken. Inspired by the electric circuit analogy method, we limit the widths of the microchannels in the 2D optimization process. To reduce the difference of fluidic performance between the 2D model and its 3D counterpart, we propose an applicable 2D optimization model, and ensure the manufacturability of the obtained layout, combinations of several morphology-mimicking filters impose maximum or minimum length scales on the solid phase or the fluidic phase. Two typical Lab-on-chip functional units, Tesla valve and fluidic channel splitter, are used to illustrate the validity of the proposed application of length scale control.


Introduction
With the trend of miniaturization in recent years, Lab-on-a-chip devices have been widely used in the area of the analysis, synthesis, and separations of biofluidics due to the advantages of high efficiency, portability, and low reagent consumption [1]. Many functions of the conventional analytical laboratory can be achieved on a centimeter-level chip, such as injection, mixing, reaction, cleaning, separation, and detection. In recent years, various microfluidic devices have been designed and applied, such as micropumps, microvalves, micromixers, and microchannels [2][3][4][5]. To pursue more functions and higher efficiency, the design and optimization of microfluidic devices is a continuously studied topic.
Borrvall et al. [6] performed topology optimization for fluidic flow governed by the Stokes equation. Their approach was later generalized to the Navier-Stokes equations by Gersborg-Hansen et al. [7] and Olesen et al. [8]. Topology optimization methods have been used to successfully design lab-on-a-chip microfluidic devices, such as microchannel splitters, micropumps, no-moving-part microvalves, and micromixers [9][10][11][12][13][14][15]. The topology optimization model of fluidic flow usually uses a simplified 2D problem. Extruded versions of the 2D optimization results typically comprise the 3D designs to be fabricated. We remark that if the 2D design is stretched so that the 3D flow path has infinite depth, then the full 3D problem describing the fluidic flow reduces to the 2D problem used in the optimization.
However, in practice, the depth of the flow channel is not infinite. The limited flow depth and the no-slip condition, which states that the flow speed is zero along the walls, together make the fluidic performance of the 3D model very different from that of the 2D model. This discrepancy significantly limits the application of the obtained 2D results. Some simplification methods have been proposed by build pseudo-3D models in recent years [16][17][18][19], and these research works are mostly focus on heat sinks problem. Zhao et al. [20] impose length scale control on the cooling channel. Based on our previous work [9,14], the length scale control methods are used to design the microfluidic devices with physical constraints (Diodicity for the example of Tesla valve, or the ratio of flowrate at the outlet for the example of splitter) in this paper.
Electric circuit analogy is a method of calculating and designing microchannel networks by analogizing the foundation parameters of the fluid to those of an electrical circuit [21]. When the height of the channel is constant, the hydraulic resistance is linear with the length and non-linear with the width of the channel. Inspired by this, we limit the widths of the microchannels in the optimization process to reduce the difference between the fluidic performance of the 2D and 3D model. Combinations of morphology mimicking filters [22][23][24] impose maximum or minimum length scales on the solid phase or the fluidic phase. The length scale control also ensures the manufacturability of the obtained layout, which avoids unreasonable small-sized islands inside the fluidic channel or thin-sized width of the fluidic channel. This paper is organized as follows-the topology optimization method and electric circuit analogy method are introduced in Section 2; the length scale control method using the morphology-mimicking filters and several numerical results are presented in Section 3; the discussion and conclusion are stated in Section 4.

Fluid Flow Model
Based on the continuity assumption, Navier-Stokes equations describe the fluid flow where f is the body forces acting on the fluid and u, p, ρ, and η are the velocity, pressure, density, and viscosity of the fluid, respectively. To model the fluid-solid interface, we let f be an artificial friction force, similarly as proposed for the Stokes flow by Borrvall et al. [6]. More precisely, the artificial friction force f = −α(γ)u, where α is the impermeability of the artificial porous material and γ is a material indicator function. The function γ varies in the interval [0, 1], where 0 and 1 denote the solid and fluid phase, respectively. The impermeability of the porous material is the interpolation function of the function γ: where α max is the impermeability of the solid phase and q is a positive value used to adjust the convexity of the interpolation function. To obtain perfect impermeability of the solid no-slip boundary, α max should be infinite; but a finite number has to be chosen to ensure numerical stability. For the topology optimization of fluid, the aim of the optimization is usually to minimize where (u) = ∇u + ∇u T /2, | * | is L 2 norm. That is, the typical aim is to minimize the viscous dissipation inside the computational domain Ω.

Electric Circuit Analogy Method
The electric circuit analogy method provides a fast calculation of laminar flow in microchannels based on the analogous behavior of hydraulic and electric circuits with correlations of pressure to voltage, volumetric flowrate to current, and hydraulic resistance to electric resistance [21]. The hydraulic resistance is an important parameter in the design of microfluidic networks by electric circuit analogy. The hydraulic resistance R H of a rectangular microchannel section can be calculated as [25] R H = 12ηL where η is the viscosity of the fluid and w, h, and L are the width, height, and length of the microchannel section, respectively. Figure 1 illustrates the hydraulic resistance changes with the length (red), height (blue), and width (green) of the channel, where the viscosity of the fluid η is normalized as 1.
When the height of the channel is constant, the hydraulic resistance changes linearly with the length and nonlinearly with the width of the channel. By limiting the width of the flow channel, the hydraulic resistance changes linearly with the length. The hydraulic resistance ratio between different parts of the fluid network can remain unchanged in different stretching heights, so the difference of fluidic performance between the 2D model and its 3D counterpart is small. In order to preserve the flexibility of topology optimization, the changes of the channel width is constrained within a certain range.

Numerical Experiments
To demonstrate the validity of the proposed method, which can reduce the difference of fluidic performance between the 2D model and its 3D counterpart by using the size control of the widths of the fluidic channels, two examples are used in this paper. The viscosity and density of the fluid are 10 −3 Pa·s and 998 kg·m −3 , respectively. The inlet velocity is chosen to have a parabolic profile. The intermediate fluid velocity in the height direction at the entrance of the 3D model is selected to be consistent with the 2D model's inlet velocity. For all numerical experiments, Svanberg's method of moving asymptotes [26] solves the resulting optimization problem.
We discretize our computational domain Ω by N quadrilateral elements. Comsol multiphysics solves the governing Navier-Stokes equation using bi-linear tensor product elements together with GLS (Galerkin least squares) streamline diffusion with the crosswind diffusion parameter set to 0.1.
We let M denote the number of nodal points in the discretization of the design domain Ω D . We define our material indicator function γ to be an nodal-wise bi-linear function whose restriction to Ω \ Ω D satisfies γ | Ω\Ω D ≡ 1. The M nodal values of the material indicator function in Ω D are obtained by using a harmonic mean based open-close f W-mean filter. For a comprehensive review of mathematical morphology, we refer the reader to the Heijmans article [27]. Figure 2 illustrates the four-basic morphological operators: dilation, erosion, opening, and closing. In colloquial terms, given a brush with shape B, then the opening of S by B holds all points the can be filled by using this brush without touching any point outside S; and given a fully filled space and an eraser with shape B, the closing of S is the region that remains after having erased as much as possible without removing any point in S. Or, in other words, applying the opening operator to a set S removes all small parts from S; similarly, applying the closing operator to a set S removes all small internals holes. The region obtained by the open operation can be written as a union of translates of B, and the complement of the region obtained by the close operation can be written as a union of translates of B. It can be shown that The use of a cascade of filters to approximate the open-close operator was suggested by Sigmund [28], and harmonic mean based filters were introduced by Svanberg et al. [29]. Using the notation for f W-filters [30], we define the equally weighted discrete harmonic erode and dilate operators with radius r > 0 and parameter β > 0 by , respectively; the neighborhood matrix G r has entries g ij = 1 if the distance between the centroids of elements i and j is smaller than r, else g ij = 0; and D r = diag{G r 1}, where 1 is the M × 1 vector with all entries set to one. The vector γ that holds the nodal values of the material indicator function γ in Ω D is where ξ is our design vector. In the numerical experiments, we use parameter β = 10 −4 . The values for r o and r c dictate the desired minimum length scale for the fluidic and solid phase, respectively. These values are selected separately for each experiment. We remark that the existence of a minimum size scale for both phases cannot be guaranteed when using open-close filters [31]. However, in many cases the optimized designs, obtained by using such filters, possess a minimum size scale on both phases [24,28].
To impose a constraint on the total volume of the fluidic phase as well as to constrain the maximum length scale of the fluidic phase, we add additional two constraints [22,23] where V * is the maximum fraction of Ω to be filled with the fluid, ε = 10 −4 , parameter r m specifies the desired maximum length scale of the fluidic channels, and v is a vector with entries where ϕ i is basis function corresponding to the ith node in Ω D and |Ω D | is the area of Ω D .

Tesla Valve
The Tesla valve is a basic device for controlling fluid motion [14,32]. Figure 3 shows the computational domain Ω for the Tesla valve. We denote the channels connected to the inlet and outlet by Ω C , and the Tesla valve occupies the domain The domain for the Tesla valve is a two-port network, where boundaries Γ in and Γ out represent the two ports. We consider the flow in the two directions from Γ in to Γ out and vice-versa. Henceforth, we denote the two flow directions by the forward and backward direction, respectively. In both the forward and backward direction, the flow profile is specified at the port from which the fluid flows. We let p f and u f denote the pressure and velocity, respectively, for the forward-directed fluid field. Similarly, p r and u r denote the pressure and velocity, respectively, for the backward-directed fluid field. The forward flow can be modeled by Navier-Stokes Equation (1) in Ω together with the boundary condition where u 0 is the given flow velocity profile at Γ in . The backward flow can be modeled by Navier-Stokes Equation (1) in Ω together with the boundary condition u r = u 0 on Γ out , [−p r I + 2η (u r )] · n = 0 on Γ in , where u 0 is the given flow velocity profile at Γ out . The performance of Tesla valve is measured by diodicity, which is defined as the ratio of pressure drop of the forward flow to that of the reverse flow. From the kinetic energy viewpoint, the diodicity is also expressed as the ratio of the consumption of the forward pressure work to the reverse pressure work, under the assumption that the flux in the two opposite directions is equal. The diodicity of the Tesla valve can be alternatively calculated as [13,14,32,33] The topology optimization problem is For the numerical simulations, each square space with an edge size equal to L is discretized by 20 × 20 quadrilateral elements. Hence 120 × 120 elements discretize the design domain Ω D . The values of α max and q in the topology method are chosen to be 5 × 10 8 and 1, respectively. Here, we set length scale parameters r c and r o small enough to ensure that the neighborhood of each node in Ω D only consists of the node itself. This yields that γ = ξ. The initial value of the design variable ξ is 1, the Reynolds number is 100, and V * is 0.55. Figure 4 shows two optimized Tesla valves. Due to physical problems and geometric constraints, there is a lower limit for the diodicity of the Tesla valve in a particular design domain. And the optimization results depend on the optimization process. Initially, we set a low C max and then increased during the optimization. The left and right valves are obtained for the value of C max is 0.351 and 0.551, respectively. We create body-fitted 2D models that only consider the fluidic phase. The diodicities computed by using these body-fitted 2D models are 0.408 and 0.600 for the left and right design in Figure 4, respectively. Next, we add a constraint on the maximum length scale in the optimization model and consider the problem min The maximum and minimum length of the fluidic phase are r m = 200 µm and r o = 100 µm, respectively, and the minimum length of the solid phase is r c = 60 µm. For this test, we set V * = 0.25 and C max = 0.5. Figure 5 shows the optimized Tesla valve with the value of C max = 0.829. We can see that the optimized design is very similar to the two commonly used Tesla valves in series. We create a body-fitted 2D model that only considers the fluidic phase. The diodicity computed by using the body fitted 2D model is 0.882. It is clear that the design constraint of diodicity is violated because of using size control during the optimization procedure.
where the subscripts 2D and 3D represent the 2D model and its 3D counterpart, respectively. In the ideal case, δ 1 is equal to 0. Figure 9 shows the difference of three optimization results between 2D and 3D models with the ratio of depth to width of the inlet n. The difference in fluidic performance between the optimized 2D model and its 3D counterpart gradually decreases as the height of the 3D model increases. This confirms that the 2D model assumption that the height is infinite, the higher the stretching height of the 3D model, the smaller the difference between the 2D model and the 3D model. The difference of the optimization results using dimension control is significantly smaller than that of unused optimization results, and it is acceptable when n is greater than 1.   n Figure 9. The difference δ 1 of three optimization results between 2D and 3D models with the ratio of depth to width of the inlet n.

Microfluidic Splitters with Equivalent Outlet Flowrate
As our second test problem, we study the problem of designing microfluidic splitters originally proposed by Zhou et al. [9]. Figure 10 shows the computational domain Ω = Ω C ∪ Ω D for this test case. The width of inlet L is 100 µm. Here, we have one inlet Γ in and three outlets Γ we denote the union of the outlets by Γ out . Given a specified velocity profile, given by u 0 on the inlet, the flow is modeled by Navier-Stokes Equation (1) in Ω together with boundary condition (11). In this example, we consider two fluidic cases in a four-port splitter. The flow direction is from Γ in to Γ out in both cases. In contrast to the previous example, we will require a specified flowrate, at selected ports of the domain. We let p real and u real denote the pressure and velocity, respectively, for the so-called (in this work) real case. The real case is modeled by Navier-Stokes Equation (1) in Ω together with the boundary condition Here Q in is the given flowrate at Γ in . Similarly, we let p ref and u ref denote the pressure and velocity, respectively, for the so-called reference case. The reference case is modeled by Navier-Stokes Equation (1) in Ω together with the boundary condition The topology optimization problem is The first term of the objective function is a fraction, in which the numerator is the sum of the viscous dissipations of the two cases and the denominator is the volume of the solid phase, to obtain an appropriate volume automatically. The second term of the objective function is the difference between the viscous dissipations of the two cases, and θ is the scalar number used to penalize this difference. By this method, the flowrate constraint of outlet in the real model is realized. The values of α max and q in the topology method are 5 × 10 9 and 1, respectively. The initial value of the design variable ξ is 0.5. During the optimization, θ increases from 0 to 100. The left column of Figure 11 shows designs optimized without length scale control. That is, just as for the results in Figure 4, r c and r o are small enough to ensure γ = ξ. The Reynolds numbers for the simulations were 0.1 (top row), and 10 (bottom row), and the value of b in the objective function was 0.1.
The right column of Figure 11 shows designs optimized under a maximum length scale constraint for the fluid domain. That is, we solve optimization problem (20) appended with constraint (9).
The maximum length of the fluidic phase is 100 µm, and the value of b used in the objective function is 0.1 and 0.2 for Reynolds numbers 0.1 and 10, respectively.
Having access to the 3D models, we compute the difference where Q (i) is the flowrate at port Γ (i) out . In the ideal case, δ 2 is equal to 0. Figure 12 shows the difference for 3D models corresponding to optimized designs in Figure 11 as functions of the ratio of depth to width of the inlet n. The difference in fluidic performance between the optimized 2D model and its 3D counterpart gradually decreases as the height of the 3D model increases. The difference of the optimization results using dimension control is significantly smaller than that of unused optimization results when Re = 10. and it is not significantly smaller than that of unused optimization results when Re = 0.1. The value of Reynolds numbers maybe influences the effect of the proposed method.

Discussion and Conclusions
We can see that the difference in fluidic performance between the optimized 2D model and its 3D counterpart gradually decreases as the height of 3D model increases in Figures 9 and 12. This conforms to the 2D model assumption is that the height is infinite.
There is a significant difference in fluidic performance in the case where no length scale control is applied. The difference in the optimization result of length scale control applied is much smaller. Imposing size control of the widths of the fluidic channels can effectively reduce the difference of fluidic performance between the 2D model and the 3D model obtained by stretching the 2D model, so that the 2D optimization results of flow have more practical significance. Moreover, the length scale control also ensures the manufacturability of the obtained. Numerical results demonstrate the validity of the proposed method. We hope that this paper will spur further interest and development on this topic.
Author Contributions: Y.G., H.P. and Z.L. wrote the topology optimization code; E.W. and Z.L. wrote the size control code; Y.G. and Z.L. performed the numerical examples; Y.G., E.W. and Z.L. wrote the paper. All authors have read and agreed to the published version of the manuscript.