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.


Introduction
Passive valves without any moving components are an intriguing technological and economical alternative to active valves. Such valves are also called Tesla valves in honor of their inventor Nikola Tesla 1 and are still the subject of current research today [2][3][4] . Two obvious advantages are the high robustness and the simple manufacturing. Even more important for applications in micro-, nano-, or even pico-fluidics is that it is intrinsically much more difficult to miniaturize valves with active components than passive designs. Unfortunately, established designs for passive valves can not simply be scaled down because they tend to loose their functionality on smaller scales with naturally smaller Reynolds numbers 5 . Miniaturization of passive valves, thus, calls for new valve designs.
The rectifying effect of Tesla valves is not based on any moving components, but relies solely on inertial forces and the resulting viscous losses, which are asymmetric for the two flow directions. This asymmetry leads to the effect that the fluid can flow rather easily through the valve in one direction, namely the forward direction, whereas the fluid experiences a significantly greater flow resistance in the reverse direction, hereinafter referred to as the backward direction. The ratio of these flow resistances describes the efficiency of a Tesla valve and is called the diodicity Di.
The development of efficient valves is an area of active research with many recent publications [5][6][7][8][9] . For practical applications, diodicities of order 2 have been reported for the Reynolds number range of Re ≈ 100. This work focuses on flow rates smaller than 20 μl s −1 , which correspond to a Reynolds number of less than or about 35, with the exact value depending on the choice of the characteristic length and velocity. This parameter range presents a particular challenge for the design of efficient passive microvalves: While for Re ≳ 100, various simulations provide sufficiently large Di ≳ 2 7-9 , simulations for Re ≈ 50 consistently find only marginal diodicities Di ≳ 1.0 5,6,10 . A paradigmatic application of Tesla valves is the rectification of the volume flow in, e.g., membranedriven 11,12 or piston-driven micropumps 13 . The idea of passive rectification can be taken a step further by combining Tesla valves with a non-mechanical pumping principle 14,15 . For the concept of such an electrokinetic pump, this type of valve is an excellent candidate, as the same manufacturing processes are used for the fabrication of the pump and the valves. While the present manuscript focuses on Tesla valves for micropumps with flows at small Reynolds numbers (Re < 35), our design strategies and objective functions can also be applied and adapted to other passive non-mechanical microfluidic devices such as micromixers [16][17][18][19] or gas decompression, e.g. in hydrogen fuel cells 20 . All these and similar applications will profit from the increased robustness and reliability, as well as the smaller footprint and inexpensive manufacturing resulting for nonmechanical static designs.
Fundamental research on the working principle of the valves was done by Bardell 21 . Gamboa et al. 17 and Truong et al. 22 worked on the optimization of Tesla valves starting from the original design of Nikola Tesla. Lin et al. 23 proposed a topological optimization method to increase the diodicity even further. Their valve designs achieve very high diodicities at medium to high flow rates. A remaining major challenge lies in the fact that the achievable diodicity of Tesla valves depends strongly on the flow rates and is usually very low for small flow rates. Tesla diodes can work in practice only when the Reynolds number is larger than one, i.e. inertia plays a role. If the Reynolds number becomes close to unity, the flow is typically reversible and the Tesla diode shows the same pressure drop in both directions 24 . In a simplified geometry (e.g. nozzle-diffusor valves) the diodicity is expected to increase linearly with Reynolds number. However, for more complex structures optimal Reynolds number ranges can be determined. Typically the massflow rate or volume flow rate is controlled in microfluidics. Both correspond to each other for liquids as in typical pressure ranges for microfluidics the liquids can be considered to be incompressible. Assuming a constant volume flow rate, channel height, and fluid properties, the Reynolds number is inversely dependent on the channel width. In Tesla valves the channel width changes with streamwise direction. Therefore, the Reynolds number should be determined at the minimum crosssection or the inlet of the Tesla diode. We decided to use the inlet as reference as also done in the literature 10,23 to be able to compare the different geometries here. Unfortunately, in microfluidic systems flow velocities are typically small, making the use of Tesla valves inefficient up to now. This problem was already recognized by Deng et al. 25 . They used a topological optimization method for the design of Tesla valves, which should be efficient for small Reynolds numbers. However, the achievable diodicities are not particularly high (of order ≈ 1.25) even with a twelve-stage valve design. Furthermore, the early research indicated the necessity of a relatively large size for efficient valves, which counteracts the purpose of microfluidics. Another approach was presented by Tao et al. 26 , who used asymmetric converging-diverging microchannels to fabricate efficient valves for small flow rates. This approach already yields high diodicities of maximum 1.77 at moderate flow rates but is still not efficient enough at very low flow rates. A comprehensive study of existing passive valves at low Reynolds numbers was performed by Fadl et al. 27 , and none of the structures studied showed a sufficiently high diodicity. In this paper, an improved multi-stage optimization procedure is presented, which generates Tesla valves that are efficient even at very low flow rates or Reynolds numbers. At the same time, these valves have a small footprint and are easy to manufacture. Different valve designs are shown and characterized by fluidic measurements. Besides the direct measurement of the diodicities, the flow profiles are determined by micro-particle image velocimetry (μPIV) measurements and compared with the simulated predictions. We start with a description of the topological optimization. While the optimization is done quasitwo-dimensional, taking into account the finite height only approximately, the subsequent simulations predicting the properties of the optimized designs are done fully three-dimensionally. Thereafter, we present a discussion and experimental verifications for both diodicities and velocity profiles. Finally, we summarize the materials and methods used.

Topological optimization
This section explains the proposed algorithm for the optimization of Tesla valves. Topological optimization 28 methods have already been proven to be useful for dimensioning efficient passive valves. A seminal paper on this subject was written by Borrvall et al. 29 . The method was subsequently refined by Lin et al. 23 . However, the method described therein does not provide efficient valve designs at small flow rates. This can be explained by the fact that for small flow rates, the pressure difference between the forward and backward directions becomes smaller, making the optimization problem more difficult. To address these issues, a multistage optimization procedure with a customized objective function was developed. In addition, we found the choice of the form of the optimization region to play a decisive role for the result. A typical variant for the possible choice of the optimization area which leads to small and efficient valve designs at the same time is presented in Fig. 1a. For the optimization, it is necessary to define what an 'optimal valve design' actually means. In the past, the requirement was primarily to achieve a high diodicity, defined as the ratio of the pressure differences between the forward and backward directions at a constant flow rate: However, valve designs that show only high diodicity are usually not applicable in a realistic microfluidic system. In addition to high diodicity, usually a low overall fluidic resistance in the forward direction is required. Further conditions such as the minimum structure size matter for manufacturing. For actual applications, the filling of the system with a liquid should be possible without forming any air bubbles inside the structures. All these and similar requirements can be accounted for by a target-oriented choice of the objective function combined with a multistage optimization procedure. For the optimization, a region is first defined where the geometry can be changed during the optimization. This area is called the design domain. Having typical microfluidic systems in mind and in view of the numerical effort, we limit ourselves here to two-dimensional design domains. In the design domain, a function γ C is introduced which describes the material distribution: A value of γ C = 0 shows that there is material in the respective point and a value of γ C = 1 means that the fluid can flow freely. In order to be able to apply a derivative-based optimization procedure, the permissible range of γ C is extended to the whole range 0 ≤ γ C ≤ 1. Nevertheless, the goal of the final design is for the function γ C to converge to either 0 or 1 at each point. The optimization task is now to find the optimal value of the γ C on each mesh element within the design domain. However, the direct use of this approach would lead to the material  distribution being strongly dependent on the size of the mesh elements used. In addition, very small and highly porous structures would possibly be obtained, which cannot be manufactured. Therefore, the material distribution γ C is filtered and, thus, smoothened by solving the Helmholtz PDE 30-32 where γ F is the filtered material distribution and R min is a minimum length scale. Good results are obtained if the filter length is chosen equal to the minimum mesh length m sz . However, this filtering leads to the fact that γ F can take values greater than 1 or less than 0. For this reason, the function is projected back to the range from 0 to 1 using a hyperbolic tangent function by solving the equation where γ P is the back-projected material distribution, γ β is the projection point and β is the projection slope that controls the projection steepness. The profile γ P is now used as input for the fluiddynamic calculations. To make the flow dependent on the material distribution, a Darcy force is introduced, where α depends nonlinearly on the projected material distribution γ P as follows: The quantity q is a freely chooseable factor that controls the convexity of the function, which must, however, be chosen appropriately in order to obtain meaningful designs. The factor α max controls how much the liquid can penetrate into the solid region. A very high value for α max is desirable to imitate the behavior of an impenetrable solid material, however a too large value of α max causes the solver to converge to a local minimum early. In the method presented here, the value of α max is adjusted during the optimization process. To be specific, α max is chosen dependent on the size of the mesh elements as where η is the viscosity of the liquid and α 0 is a freely chooseable parameter. Borrvall and Petersson 29 have shown that direct maximization of the diodicity according to Eq. (1) constitutes an ill-posed optimization problem. Instead, they proposed to optimize the ratio of the energy dissipations for the two flow directions where the Einstein summation convention is used, u i are the components of the velocity field, and the indices B and F refer to the backward and the forward direction, respectively. Lin et al. 23 showed under which assumptions both objective functions lead to the same result. The Borrvall-Petersson approach is also used here, but the objective function is extended by further terms. The flow profiles entering, e.g., Eq. (6) are obtained by solving the laminar Navier-Stokes equation twice in each optimization iteration. The term Àαðγ P Þ u ! corresponds to the Darcy force, which couples the material distribution to the flow profile. The last term À12ðη=d 2 z Þ u ! results from a so-called shallow-channel approximation 34 . This term is introduced in order to at least partially take into account the influence of the low channel height during the two-dimensional simulations used for the optimization. For small flow rates, our experience shows that the direct use of Eq. (6) as an objective function is not effective. Therefore, the following adjusted objective function is used: This objective function now consists of three different terms. The parameters w 1,2,3 ∈ [0, 1] are weighting factors for the different, partly competing design goals. The first term corresponds to the ratio of the energy dissipation according to Eq. (6). The second term corresponds to the actual optimization of the diodicity according to Eq. (1).
Since the optimization of energy dissipation does not exactly correspond to the optimization of diodicity, this term is helpful to guide the optimizer in the right direction. The third term favours designs that fill as much of the structure as possible with material in order to ease fabrication: A O is the total area of the design domain.
Together with the filtering of the material distribution via Eq. (3), this term prevents structures that are too small for actual fabrication.
Using this objective function, a multi-stage optimization procedure is carried out, which adjusts the value of the minimum mesh length m sz in each step. First, the value of m sz is chosen rather large and then gradually reduced. This leads to the effect that there are fewer mesh elements at the beginning and the search space is thus comparatively small. In addition, the liquid can according to Eq. (5) penetrate the material stronger, which increases the robustness of the optimization at the beginning. In each subsequent optimization step, the previously found material distribution is used as a starting point for further optimization and the mesh size is reduced. A possible problem with this method is that only a single flow rate _ V 0 is taken into account during optimization. Thus, a single operating point of the valve must already be determined at the beginning. This leads to the fact that the diodicity can be very high for the specific flow rate but at the same time, the response curve is usually steep and drops off quickly for smaller flow rates. However, in many cases, valve designs that work efficiently over a wide range of flow rates are desirable. In addition, a strongly non-linear response curve leads to a strongly non-linear system behavior, which is often not wanted. For these reasons, a further optimization step with an adjusted objective function is carried out within our multi-stage optimization. This last stage uses the following objective function where the optimization is now carried out over a given range of flow rates simultaneously: Here, the set of flow rates _ V i , i = 1, …N is specified in a given range of interest. Since an efficient valve design is already available at this point of the multistage optimization, the diodicity instead of the ratio of the energy dissipation can now be maximized directly, which is expressed by the first term. In order to increase the influence of the usually small diodicities at small flow rates to the objective, the maximum flow rate _ V max divided by the i-th flow rate _ V i is used as a weighting factor for the diodicity. In addition, another term is introduced which favours the minimization of the pressure difference in the forward direction. In contrast to solely optimizing the diodicity, not only the ratio of the pressure differences but also a small total fluidic resistance plays a relevant role. Deng et al. 25 showed the usefulness of this additional optimization term. To weight these terms similarly, the pressure difference is divided by the inlet velocity u in for the i-th flow rate multiplied by the density of the liquid ϱ. In addition to the choice of objective functions, the form of the design domain plays a decisive role for the result. The form of the design domain used here is shown in Fig. 1a. This curved shape of the design domain was chosen because the normally used sharp-edged design domains, see e.g. Lin et al. 23 , lead to air inclusions and bubbles during the filling process of actual manufactured structures. In summary, this study presents the following four novelties in order to find designs of efficient microfluidic Tesla valves at low flow rates that are still easy to manufacture: 1. A customized objective function that simultaneously maximizes the diodicity and energy dissipation and also prevents the formation of small structures. 2. A multi-stage optimization procedure in which the element size of the underlying mesh is reduced successively. 3. A final optimization step in which the diodicity is simultaneously maximized over a certain flow rate range and, if desired, the fluidic resistance is also minimized in the direction of flow. 4. An appropriate choice of the shape of the design domain, which favors designs that allow easy manufacturing and a bubble-free filling process.

Simulation results
In this section, exemplary simulation results and the parameters used are presented. All simulations were carried out using the finite element software COMSOL Multiphysics ® (the Supplementary Information provides more details about the implementation in COMSOL). For the optimization, the gradient-based Sparse Nonlinear OPTimizer SNOPT 35 is used. The meshing is done with triangular elements. Constant shape functions are used to approximate the material distribution function γ P and first order elements are used to discretize the velocity and the pressure field. At the beginning, the material distribution on the entire design domain is initialized with a value of 1, which means that there is no material within the design domain at the beginning. A projection slope β = 8 and the projection point γ β = 0.5 are used. A value of q = 1 is set for the Darcy penalization. In the first step, a certain number of optimization steps are carried out while using a single fixed flow rate value. Each optimization step consists of a given number of optimization iterations. Our experience suggests that a four-stage optimization procedure is a good trade-off between computational effort and valve efficiency. Between two successive optimization steps, the size of the mesh elements is reduced step-wise by 5 μm from an initial value of 20 μm to 5 μm. A maximum of 500 iterations are performed in the first four optimization stages. Following the four-stage optimization using the objective function defined by Eq. (8), a further optimization step is carried out using Eq. (9) as objective function. A maximum of 1000 iterations are performed in this last optimization step. For all simulations presented in this manuscript, α 0 was set equal to 25. Trying different values, we found that good results are obtained for w 1 = w 2 = 1 and w 3 = 0.05. Examples of the simulation results after each step of the optimization are shown in Fig. 2a-f (as Supplementary Information a video is available, which shows the material distribution during the optimization). It can be seen that a diode-like structure is formed at an early stage. In the progress of the optimization, additional island structures already known from conventional topological optimizations are added to the structure. One notices that many small additional structures are formed by the software before the final optimization step. However, these are removed again automatically in the last optimization step, as they do not improve the diodicity at small flow rates. It becomes clear that the optimization process favors nozzle-like structures. These nozzles cause the flow velocity to be greatly increased locally. Together with the 'cleverly' placed obstacles, this greatly increases the flow resistance in the backward direction.
Three different valve structures for different parameter combinations are shown in Fig. 1b. The parameters used for the optimization are summarized in Table 1. For the optimization of design V A and V B , all weights except w 4 were chosen identically. The resulting performance difference highlights the influence of the additional minimization of the pressure difference in the forward flow direction in the final optimization step: Both designs show  Figure f shows the result after the final optimization step using the objective function of Eq. (9). A three-dimensional model can then be obtained from such a design by extrusion, whereas the two-dimensional geometry is sufficient for manufacturing by means of lithographic processes.
only minimal geometric differences leading to an also almost identical diodicity. However, the pressure difference for the two flow directions is significantly lower for design V B as shown in Fig. 1c.
The described optimization algorithm yields twodimensional geometries at first. Real manufactured valve structures, on the other hand, are always threedimensional and, in the case of microfluidic systems, the aspect ratio of channel width to channel height is usually very high. Even though the additional term À12 η Eq. (7) already partially accounts for the finite channel height during optimization, an exact agreement in the three-dimensional case with the optimization results cannot be expected. In order to simulate the diodicity as accurately as possible, the two-dimensional material distributions were translated into three-dimensional geometries. For this purpose, the isoline γ P = γ th is determined where γ th is a threshold value. The resulting curve provides closed areas which are then extruded and subtracted from the extruded initial simulation area. In all results shown below, a value of γ th = 0.5 was chosen. The threshold value plays only a minor role for the final geometry, since the presented optimization algorithm provides very sharp transitions in the material distribution. Finally, for the three-dimensional valve geometries, the flow ratedependent diodicity can be calculated numerically: For this purpose, the Navier-Stokes equations ϱð u ! Á ∇Þ u ! ¼ À∇p þ η∇ 2 u ! are solved for both flow directions and the pressure differences are determined. All the results shown subsequently, as well as the pressure curves in Fig. 1c, were determined for the three-dimensional geometries. To validate the simulation results, the dependence of the shape of the isolines on the choice of the threshold value γ th and the dependence of the results of the threedimensional simulations on the mesh size were investigated. The results of the convergence study are summarized in the Supplementary Information.

Pressure measurements and diodicity
In this section, measurement results gained from fabricated test chips are compared with simulation results. The diodicities were determined by measuring the pressure drops for the both flow directions as described in the methods section of the diodicity measurement. An example for the pressure drops for a range of flow rates is show in Fig. 5b. As expected, the flow resistance increases with an increasing flow rate. The typical non-linear behavior of the Tesla valves is apparent. In addition, it can be seen that the backward direction shows a significantly higher flow resistance. The measured pressure curves can be used to directly calculate the diodicity. However, since the pressure was not measured directly at the outlet or inlet of the valves, there were additional pressure drops Δp ind independent of the flow direction. In this case, the measurable diodicity is given by: In the simulations, the contribution Δp ind is zero by definition and Eq. (10) reduces to Eq. (1). Since the term Δp ind in the measurement is small but not zero, an exact agreement with the simulation results cannot be expected, and indeed the measured diodicity turned out to be smaller than the simulated diodicity. Figure 3a shows the comparison between the measured and simulated diodicities for the different valve designs. For a given geometry the Reynolds number depends only on the volume flow rate and the upper x-axis presents the inlet Reynolds number. Even though the measured diodicities are, as discussed, usually slightly smaller, there is very convincing agreement between the simulation and the measurement results. It can be seen that a minimal diodicity of 1.3 can be reached for all geometries for inlet Reynolds numbers between 10 and 35 for three different valve designs V A , V B , and V C . Valve designs V A and V B show a rapid increase in diodicity for small flow rates, which is subsequently followed by a plateau of the diodicity. In contrast, the diodicity for design V C increases slower with Reynolds number (or volume flow rate) but finally reaches significantly larger values. This can be explained by the fact that the final optimization step used a larger flow rate range (see Table 1). The results show that a measured diodicity of approximately 1.8 is achieved at a flow rate of only 20 μl s −1 (corresponding to a Reynolds number of approximately 36 calculated at the inlet), which to our knowledge is currently the highest diodicity at such low flow rates and ΔDi = Di − 1 is between two to three times as high as the diodicities reported before. At the same time, the absolute pressure drop is low, making the valves ideal for microfluidic applications.

Velocity field measurements
In addition to the measurement of the pressure differences for the determination of the diodicities, the velocity fields were measured using two-dimensional μPIV with the setup described in the methods section of the velocity measurement. The measured velocity fields can be directly compared with the simulated velocity fields and they show a good qualitative agreement. In Fig. 2b, c, the measured velocity is shown for design V A in the upper row of the subfigures for forward (b) and backward (c) flow. The data for the other designs are presented in the Supplementary Information. What becomes obvious already from the visual inspection is that the magnitude of the velocity in the experiment is always a bit higher than in the simulation. The mean ratio between the measured and simulated velocities is 1.45. The higher velocity can be explained due to the actual channels being smaller in height compared to the numerical simulation. The smaller height is caused by the etching and the processing of the Ordyl bond medium that typically shrinks during the bonding or becomes thinner during the development process. As the volume flow rate was controlled in the experiment a smaller cross-section results in larger velocities. To quantitatively assess the agreement in the patterns a correlation between the experimental and simulated velocity fields was performed (for more details see Sachs et al. 36 ). This analysis gives correlation coefficients 37 of ≥95% for the data of Fig. 3 for both directions, where a correlation of 100% would stand for a perfect agreement. As the flow topology in the investigated Reynolds number range does not change significantly the μPIV-experiments validate the numerical simulations and the whole optimization approach can be expected to result in optimized geometries even if small tolerances in the final velocity field due to the fabrication arise.
In the paper, a robust topological optimization procedure was introduced for designing Tesla valves which are efficient even at low flow rates. To the best of our knowledge, the proposed process provides the most efficient valve designs currently available in this flow-rate range. The resulting designs show only large free-standing island structures and are thus comparatively easy to manufacture and simultaneously offer a small footprint. Thus, they are promising candidates for passive valve designs in microfluidic applications with low flow rates. The measurements of the diodicities, as well as the flow profiles on three-dimensional test structures, show good agreement with the simulation results, which confirms the predictive accuracy of the simulation. Since objective functions with different terms were used for the optimization, the final design can be influenced by the appropriate choice of weighting parameters. It should be noted explicitly that the presented multi-stage optimization procedure can also be applied for designing Tesla valves which are efficient at higher flow rates and at the same time simple to manufacture. We also stress again that only readily available and economically viable fabrication steps were used for the valves produced for experimental verification.

Fabrication of Tesla valves
The Tesla valves used for the experimental characterization and verification of the optimized designs were created in a silicon substrate using a deep reactive ion (DRIE) etch process 38 . The individual steps of the manufacturing process are shown in Fig. 4a. Due to the very high aspect ratio requirement needed to create the island structures of the Tesla valves, a specially optimized etching process had to be developed. The optimization employed statistical design of experiment (DOE) 39,40 . The DRIE process is conventionally known as a two-cycle process consisting of an etch step and a passivation step, where the etch step is based on a DC Bias to physically etch the passivation layer on the bottom of the structures by accelerated ions. In contrast, the DRIE process used here consists of three steps namely, a passivation step, a breakthrough step, and a purely chemical etch step. The breakthrough step carries out the etching of the passivation layer and this step was optimized with DOE whereas standard process parameters were used for the other two steps. The etching is done at a PlasmaPro 100 Estrelas. The corresponding breakthrough parameters namely the platen power and the breakthrough time and the results are summarized in Table 2. Using the optimized etching process 'Run 3', island structures with lateral dimensions as small as 10 μm can be etched to a depth of 200 μm. Examples of etched Tesla valves are shown in Fig. 4b. In addition to the etching of the Tesla valves, a second DRIE process was employed to etch through holes from the backside of the substrate to create openings for the fluidic connections.

Wafer-bonding
The Tesla valves are sealed with a lid by bonding with a glass substrate. A transparent lid is used to allow measurements with microscopic particle image velocimetry. The island structures of the Tesla valves pose a particular challenge for the bonding. For this reason, the 30 μmthick permanent dry film resist Ordyl SY330 was used as bonding material. In view of future applications, it should be noted that this material is biocompatible and features high chemical and thermal resistance and for these reasons has already been successfully used for the construction of microfluidic devices [41][42][43] . The dry film resist was laminated by a temperature of 105°C, a pressure of approximately 0.5 bar, and a lamination speed of 1 m min −1 over the Tesla structures and the channels.
The dry film resist that behaves like a negative resist was then structured using a conventional photolithographic process. The resist was exposed in a MA8 mask aligner from SüSS MICROTEC with a dose of 200 mJ cm −2 and developed in the Ordyl Developer XFB for a total duration of 6 min, whereby the resist only remains on the island structures of the valves and on the surfaces surrounding the fluid channels. The base and the lid wafers were then stacked together and bonded with a force of 6 kN and a temperature of 150°C where the resist cross-links and stabilizes. To achieve complete cross-linking, the temperature of 150°C was maintained for at least 20 min. The bonding process is visualized in the Supplementary Information. The final chip was proven to be watertight at least up to a water flow of 100 μL s −1 , and the bonding interface was durable even against solvents such as acetone and isopropanol. After bonding was complete, the wafer stack was diced into test chips, with four valve structures on each chip. In the final step, nanoports were attached to enable a fluidic connection of the chips. An example of a finished test chip is shown in Fig. 4c.

Diodicity measurements
This section describes the experimental setup for the measurement of diodicity of the Tesla valves as defined by Eq. (1). Thus, for the direct measurement of the diodicity, it is necessary to measure the inlet and outlet pressure of a Tesla valve for both flow directions and calculate the ratio of the respective pressure drops. To generate a constant flow rate, a high precision neMESYS syringe pump from CETONI was used. At the same time, the pressure was measured directly before and after the fluid connections with two precise MPS-pressure sensors from Elveflow. The experimental setup is shown in Fig. 5a. The flow rate was increased by 1 μl every seven seconds from 0 μl to 20 μl. Since the pressure was initially unstable after every flow rate variation, the measured values during the first and last second of each seven-second-long measurement interval were ignored. The mean value of the pressure values was calculated from the values of the remaining five seconds of each measuring interval. The resulting mean values are shown in Fig. 5b. After measuring the pressure drop for both directions, the diodicity was calculated using Eq. (1) (see Supplementary Information for more details).

Velocity measurements
For a detailed comparison of the theoretically designed structure and the actually obtained Tesla valves, velocity measurements were carried out using microscopic particle image velocimetry (μPIV [44][45][46]. For the potential of this method and recent advancements, the interested reader is referred to the results of the latest PIV challenge 47 . Small tracer particles (d P = 1μm) doped with a fluorescent dye were suspended in de-ionized water as the working fluid. The volume flow rate was adjusted to 10 μl s −1 using a neMESYS syringe pump. The region of interest was observed with an inverted microscope (Zeiss Axio Observer by Zeiss GmbH) with a magnification of M = 5. A double-pulse Nd:YAG laser with wavelength of 432 nm (Innolas GmbH) was used to illuminate the particles. A dichroic mirror and a long-pass filter allow only the fluorescent light of the particles to pass to the camera sensor. Each experiment recorded 1000 doubleframe images with a 12 bit Sensicam QE camera (PCO GmbH) with 1376 × 1040 pixel using a frame rate of~4 Hz. The time delay between frames was in the order of 10 μs and was adjusted to the specific setup to give a mean particle image displacement of about 15 pixel. A calibration target (Thorlabs GmbH, 100 μm spacing) was used for calibration and gave a calibration factor of 1185 pixel/mm. For time synchronization and image evaluation the software DaVis (LaVision GmbH) was used. To reduce the effect of the out-of-focus particles for the cross-correlation evaluation, image pre-processing was applied (subtraction of sliding minimum for background removal, min-max intensity and bandwidth filter) 48 . As only small particles were used and the seeding concentration was kept moderate, the sum-of-correlation approach 49 was used to give reliable results with high resolution. A multi-grid window deformation approach with decreasing interrogation window size and 50% overlap was applied to determine the large velocity range with low relative uncertainty with a final vector spacing of 6.75 μm in each direction. The experimental setup for the measurement of diodicity and exemplary pressure curves. a A syringe pump is connected to the first pressure sensor and then to the inlet of the valve. The second pressure sensor is attached to the outlet of the valve. At the end of the tubes, a small 3D-printed water sink is attached to allow the liquid to flow out without bubbles. Otherwise, due to the Laplace pressure, the bubbles would strongly influence the pressure signal. The analogue sensor signals are digitised using custom-built electronics and are automatically recorded. b Example of the measured and simulated pressure difference as a function of the flow rate for valve V C for the different flow directions.