Highly efficient passive Tesla valves for microfluidic applications

A multistage optimization method is developed yielding Tesla valves that are efficient even at low flow rates, characteristic, e.g., for almost all microfluidic systems, where passive valves have intrinsic advantages over active ones. We report on optimized structures that show a diodicity of up to 1.8 already at flow rates of 20 μl s−1 corresponding to a Reynolds number of 36. Centerpiece of the design is a topological optimization based on the finite element method. It is set-up to yield easy-to-fabricate valve structures with a small footprint that can be directly used in microfluidic systems. Our numerical two-dimensional optimization takes into account the finite height of the channel approximately by means of a so-called shallow-channel approximation. Based on the three-dimensionally extruded optimized designs, various test structures were fabricated using standard, widely available microsystem manufacturing techniques. The manufacturing process is described in detail since it can be used for the production of similar cost-effective microfluidic systems. For the experimentally fabricated chips, the efficiency of the different valve designs, i.e., the diodicity defined as the ratio of the measured pressure drops in backward and forward flow directions, respectively, is measured and compared to theoretical predictions obtained from full 3D calculations of the Tesla valves. Good agreement is found. In addition to the direct measurement of the diodicities, the flow profiles in the fabricated test structures are determined using a two-dimensional microscopic particle image velocimetry (μPIV) method. Again, a reasonable good agreement of the measured flow profiles with simulated predictions is observed.


Simulation and optimization details
Some important details of the implementation of the simulation in COMSOL Multiphysics ® are presented below. COMSOL specific terms and concepts are written in italics.

Optimization
The complete multi-stage optimization is performed using a single COMSOL file. Two laminar fluid flow interfaces are added to simulate the fluid flow for the forward and backward direction separately. A fully developed flow boundary condition is set at the inlet, where the mean flow velocity U in =V /A is specified. A pressure of 0 Pa is set at the outlet with the suppress backflow option selected. A no slip condition is set at all outer edges. To enforce symmetry, a slip condition is set at the center edge. The Darcy force is introduced as a volume force on the design domain. For each optimization step, a triangular mesh is introduced, where the size of the triangles on the design domain is limited by the value m z as the maximum element size. A study step is introduced for each optimization step and for the calculation of the initial flow. In each study step, the corresponding mesh is assigned. The initial values of variables solved for are taken from the previous study step to enforce a continuous optimization. The parameters which are introduced as local variables and the objective functions are assigned to each study step by using the modify physics option. A nonlocal coupling operator is introduced to calculate the energy dissipation inside of the design domain and the mean pressure on the inlets, respectively. For the calculation of the pressures, two average operators are assigned to the respective boundaries of the geometry which correspond to the inlets. In order to derive a geometry from the calculated material distribution function γ P , a filter data set with an upper bound of 0.5 is introduced in the results section. Thus, this data set contains a mesh representing the the solid structures. From this data set, a mesh part can be derived, which is exported as a file and serves as input for the subsequent three-dimensional fluid flow simulations.

Characterization of the valve structures using 3D-simulations
To create the three-dimensional COMSOL geometry, a two-dimensional work plane is added first. The mesh part from the optimization is imported onto this plane. If very small structures are present, they are removed manually. In addition, any existing interrupted line elements are repaired. Then the 2D geometry is extruded with a height of h/2, where h is the total height of the valve structures. The individual steps for transferring the two-dimensional geometry into a lithography mask or a 3D simulation geometry are shown in Fig. 1. Two laminar flow interfaces are again added to calculate the flow-ratedependent diodicity of the three-dimensional valve structures. The boundary conditions are chosen corresponding to the two-dimensional simulation. To save computing time, only a quarter of the actual valve structure is simulated, with the mirror symmetry in both directions implemented by slip conditions on the symmetry planes. In order to determine the diodicity for different flow rates, an auxiliary sweep is performed. The resulting geometry from step (a) is triangulated. For mask preparation, this triangulation is translated directly into a DXF-file using an in-house developed Matlab program.
To do this, the outer border of the triangulation is determined. The associated edges can be found by the fact that they are only part of a single triangle. Subsequently, the associated points are written to a DXF-file, which can be read directly by conventional CAD programs. In our case, Autocad is used for this purpose. If very small domains are present in the triangulation, they are finally removed in Autocad. (c) For the simulation of the valve performance, a three-dimensional geometry must be created. For this purpose, the mesh from step (b) is converted into a surface in COMSOL. If very small domains are present (see inset), they will be removed. The resulting geometry is shown in (d). The resulting surfaces are extruded. The three-dimensional valve structure is shown in (e). In addition, the channel is created as a 3D block. Then the valve geometry is subtracted from this block using a Boolean operation. (f) Final geometry, which is used for the simulation of valve performance. Only a quarter of the geometry is needed, since the top and the inner side represent a symmetry surface. The material distribution is described by the function γp. This function can take all values between 0 and 1. To derive a geometry from this, a threshold value γ th must be defined. The geometry is given by the surfaces for which the condition γp ≤ γ th is satisfied. Since the function is continuous, the boundary of these areas shifts slightly for a different choice of threshold value. However, the final mesh used in the topological optimization is very fine and the material distribution is well converged. Thus, the dependence of the position of the boundary on the choice of the threshold value is small and negligible compared to manufacturing tolerances. This is clearly demonstrated by the subfigures (a) to (c), where a different threshold value was used in each image: (a) γ th = 0.4; (b) γ th = 0.5; (c) γ th = 0.6

Mesh-independence study
The numerical results always show some dependence on the underlying discretization. Therefore, in order to make the results trustworthy, we performed two different convergence studies. First, we studied the geometry dependence of the two-dimensional valve geometry on the choice of the threshold value γ th . Since the material distribution is given by a continuous function, the edges are not ideally sharp. However, the mesh used in the final step of the optimization is so fine that the choice of the threshold value has only a very small influence on the final geometry. This is demonstrated in Fig. 2. In addition to the choice of threshold value, we also investigated the influence of the mesh size on the 3D simulations. Since the results of the 3D simulations are used for the final comparison between experiment and simulation, a good convergence must be achieved. The results of the convergence study are presented in Fig. 3. It is clear that the mesh used for the simulations is sufficiently accurate, and further refinement does not provide noticeable further improvement.

Measurements
Additional measurement results on the valves which were not discussed in the main text, are presented afterwards.

Diodicities
In Fig. 3a of the paper the measured values of the diodicity are shown as a function of the Reynolds number for the different valve designs. The Reynolds number was calculated via: where ϱ is the density of the liquid, d H is the hydraulic diameter, U =V /A in is the average velocity at the inlet, which is calculated as the ratio of the flow rateV and the area of the inlet A in and η is the viscosity of the liquid. The hydraulic diameter is calculated via: where h C = 130 µm denotes the etch depth of the valve and l C = 1000 µm is the width of the inlet. Thus, the Reynolds number is calculated as usual for the characterization of Tesla valves, see e.g. [1]. However, we would like to point out that the flow velocity shows a strong spatial dependence and is strongly increased within the nozzle structure. Thus, a correct definition of the Reynolds number is difficult, as already discussed by Tao et al. [2]. Each measurement point in Fig. 3a of the main paper corresponds to the average of three individual measurements. The diodicity Di and the standard error SE of the diodicity are calculated as described by Dunlap et al. [3] using: where σ P denotes the standard deviation and ∆p is the mean value of the measured pressure values at the respective flow rate. To more accurately quantify the error between the simulation and measured data, the squared bias SB, the mean squared variation M SV , the mean squared deviation M SD and the correlation coefficient r were calculated for the different valve designs. The calculation is based on the following equations [4]: where Di M denote the measured and Di S the simulated values for the diodicities. The quantities with a bar denote the respective mean values. The squared bias SB indicates whether our simulation model can describe the measurement well. The smaller the value of M SV , the better the simulation can explain the deviations of the measurement from the mean. Both quantities add up to the mean squared deviation M SD. Thus, for a good agreement between measurement and simulation, the two quantities SB and M SV must be small at the same time. In addition, an r-value closer to 1 indicates a better agreement between the measurement and the simulation. The results for all the quantities are summarized in Table 1 for the different designs. They have also been mentioned in the caption of Fig. 3. The equations were added to the supplemental. It becomes clear that both SB and M SV in each case are low and that the r value is close to 1. Thus, we can speak of a good agreement between the measured and the simulated data. However, it also becomes clear that the deviations between the simulation and the measurements slightly increase with increasing diodicity.    The resist is stiff enough to span the channels. (d) Using conventional lithography techniques, the Ordyl is exposed and subsequently wet-chemically developed. The resists remains only on the island structures as well as around the channel.

Microscopic particle imaging velocimetry
(e) A glass wafer is then bonded to the Ordyl layer at a temperature of 150 • C and a force of 6 kN.

Fabrication
To further clarify the process of bonding using the Ordyl SY330, the process is illustrated in Fig. 6. It should also be mentioned that the the process described here enables also the bonding of two structured substrates. This is particularly relevant in the manufacturing of micropumps.

Supplemental Video
A video is available online showing the optimization process for valve V B . At the beginning of each optimization stage, the used mesh is shown and then the calculated material distributions γ P ≤ 0.5 are visualized after every 10th iteration. The upper half of each frame shows the flow profile in the backward direction and the lower half shows the flow profile for the forward direction. The same color scale is used for all frames of each optimization stage.