Reconstruction of Tesla micro-valve using topological sensitivity analysis

Abstract: In this paper, we deal with topology optimization attributed to the non stationary Navier-Stokes equations. We propose an approach where we analyze the sensitivity of a shape function relating to a perturbation of the ow domain. A numerical optimization algorithm based on topological gradientmethod is built and applied to the 2D Tesla micro valve reconstruction. Some numerical results con rm the e ciency of the proposed approach.


Introduction
Tesla valves are no-moving-part valves that utilize uidic inertial forces to inhibit ow in the reverse direction. It was patented in 1920 by Nikola Tesla as a "Valvular conduit" [1] (see gure 1), and has since made the subject of various applications in micro-satellite [2], drug delivery [3], microbiology [4,5] and hydrocephalus treatment in medicine [6,7].  [8] The Tesla micro-valve performance is evaluated by the diodicity parameter (represents the ratio of the pressure drop in backward and forward direction) which evaluates the ability of allowing forward ow while inhibiting the reverse one, Di = ∆p backward ∆p forward .
Di erent works have focused on the optimal shape of the tesla micro-valve. However, the majority of works concerns stationary Partial Di erential Equations (PDE). Forster et al. [8] proved the possibility of using Tesla valves in micro-uidics and determined experimentally the diodicity for Reynolds number (Re) 180. Truong et al. in [9] derived numerically the optimum geometry of Tesla valve for < Re < with better diodicity than [8]. Bardell et al [10] analyzed the mechanism of the diodicity and proposed a Tesla valve optimal design for low Re. In the case when Re = , they nished with Di = . . Gamboa et al. [11] optimized the shape of Tesla valve for application with piezoactuated plenums. The obtained uid domain related to Re = is characterized by a diodicity number Di = . . In 2008, Pingen et al. [12] used the Lattice Bolzmann Method for the optimization of a micro Tesla valve without any information on diodicity. The used objective function was the pressure drop between inlet and outlet. After that in 2010, Lin et al. [13] used a topology optimization technique based on the power dissipation energy [14] of forward ow as objective function and the diodicity was built into the model as a constraint. For Re = , they found a new design of Tesla valve given Di = . . Next in 2015, Lin et al. [15] solved the Tesla valve topology optimization using the approach of material distribution with inverse diodicity as objective function and uid volume fraction as the constraint.
Until recently, there were no investigations dealing with the non-stationary case. We propose in this paper a new reconstruction method using the sensitivity analysis approach [16][17][18][19] for a non stationary ow.
The principal results of this work concern both theoretical and numerical aspects associated with the Tesla micro-valve problem. The theoretical part is related to the analysis of the topological sensitivity for the non stationary Navier-Stokes equations. The numerical part concerns the 2D optimization of the Tesla microvalve shape. The optimal shape is constructed by inserting obstacles in the considered initial domain. We build a simple and fast numerical reconstruction algorithm based on the topological gradient technique. The e ciency of the presented approach is con rmed by some numerical tests.
The paper is presented as following: Firstly we formulate the problem in section 2. Section 3 concerns the theoretical aspects. The numerical aspects are given in section 4. Finally section 5 includes Theorems proofs.

Problem formulation
Let Ω ⊂ R d , d = , a bounded domain with regular boundary Γ = ∂Ω. We consider the blood as an incompressible viscous uid ow described by the non stationary Navier-Stokes equations [20]. The velocity w and the pressure p satisfy the following system: where ν is the kinematic viscosity coe cient, G is the gravitational force, T is the computational time and w d is a given Dirichlet boundary data. Because of the divergence free condition on w, w d must necessarily satisfy the compatibility condition, where n is the unit outward normal vector along Γ.
The topological sensitivity method idea is to study the variation of a given shape function j relating to a perturbation in the uid ow domain geometry.
In structural shape optimization case (respectively electromagnetism and uid dynamics cases) a geometry perturbation means removing some material (respectively the insertion of an obstacle). Let Oz, ε = z + εO, a small obstacle inserted in Ω characterized by its center z, its size ε and its shape O. O is a bounded domain of R d containing the origin and ∂O (its boundary) is connected and piecewise C .
The shape function variation is written where ε → ρ(ε), a positive scalar function going to zero with ε z → δj(z), called the topological gradient, describes the shape function variation when an obstacle is inserted in z. It plays the role of descent direction in the algorithm of optimization.
To our knowledge, the majority of works leading with topological sensitivity method concern the stationary case such as Stokes problem [16,18], quasi-Stokes [19], stationary Navier Stokes problem [17]. We extend this method to the nonlinear unsteady Navier Stokes ow. To overcome the di culty due to the non linear operator and its associated adjoint problem we extend the perturbed velocity by zero in the inclusion which permits to use the adjoint method in the whole domain. For the time dependent term we will use the fundamental solution of the non stationary Stokes operator and decompose the velocity variation.

Main results
We deal in this section with the non stationary Navier-Stokes topological sensitivity relating to the domain perturbation. We consider the shape functions verifying the assumption (A).

. Asymptotic behavior of the velocity variation
We rst study the in uence on the velocity vε = wε − w of inserting a small obstacle Oz,ε in Ω. From (1) and (3), it is straightforward to show that (vε , pv ε ) satisfy the system We will distinguish in the following the D and D cases.

. . Three dimensional case
Theorem 3.1. There exists c > independent of ε, such that with U j is solution of (exterior Stokes problem) with {e j } j= , , is the R canonical basis.
We show by using a single layer potential (see [22]) that where with r = y , er = y r , e T r is the transpose of er and η j ∈ H − / (∂O) is a solution of the boundary integral equation Using Theorem 3.1 we obtain the following corollary.

Theorem 3.3. There exists c > independent on ε, verifying
represents the fundamental solution of the Stokes System in R with r = y and er = y r .
Using Theorem 3.3 it follows the velocity estimation in the perturbed uid ow domain.

. Asymptotic behavior of the shape function
The topological sensitivity analysis for the non stationary Navier-Stokes operator in three and two dimensional cases is given in this section. The presented results are satis ed by all shape functions j de ned by (2) and Jε veri es the Assumption (A).

. . Three dimensional case
Theorem 3.5. If Jε satis es the Assumption (A) with ρ(ε) = ε, then j de ned by (2) veri es where the matrix M O is given by u is the solution to the adjoint problem Corollary 3.6. If O = B( , ) (the unit ball), η j (y) = − ν e j , ∀y ∈ ∂O and . . Two dimensional case Theorem 3.7. If Jε satis es the Assumption (A) then j de ned by (2) veri es where u is the adjoint state solution to the problem (9).
The proofs of Theorems 3.1, 3.3, 3.5 and 3.7 are relegated to section 5. The variation δJ depends on the shape functions expressions. Some useful examples in numerical applications will be presented in section 3.3.

. Shape function examples . . First example
We de ne the shape function where W d ∈ L ( , T; H (Ω)) is a datum representing a desired uid ow state. This example concerns the L -norm shape function that has been used in geometric control problems like the optimization of location of some obstacle in a tank to approximate an object ow W d (see [16]).

. . Second example
We de ne the shape function which corresponds to the dissipation energy minimization where W d ∈ L ( , T; H (Ω)) is a given datum. It was used in several optimization problems such as minimum drag problem [23], pipe bend design [10? ], cavity example [24], reconstruction of Tesla valve [13].

Numerical results
In this section, we deal with some numerical applications to validate the obtained theoretical results given in section 3.

. Validation of the asymptotic expansion
To establish the numerical validation of Theorem 3.7, we consider the variation relating to ε of We expect to prove numerically that ∆z(ε) satis es the previously derived theoretical estimate To this aim, we consider the following data: of the considered obstacles are arbitrary chosen (see Table 1). − The shape function j is de ned by the semi-norm where W d is a given velocity state.
In this case, the function ∆z i (ε) is de ned by (see Theorem 3.7 and Proposition 3.9): The validation algorithm uses the following steps: The validation algorithm:
It corresponds to the slope of the line approximating the variation ε → log(|∆z i (ε)|) relating to log − log(ε) for each obstacle z i ε , i = , .., .
From the plotted curves in Figure 2, one deduce the slopes β i , i = , ..., in table 2.
We deduce that the numerical results con rm the behavior predicted by the theoretical estimate ).

. The Tesla micro-valve application
The hydrocephalus treatment is a very important application in medicine. The problem is to optimize numerically the design of the 2D Tesla micro-valve at Re=100. To solve this problem we consider the objective function as the forward energy dissipation and the diodicity as a constraint. The optimal domain is constructed through the insertion of some obstacles in the initial one. The problem leads to optimize the location of obstacles.

. . Shape optimization problem
We de ne Ω as the pentagon [15] having one inclined inlet Γ in and one horizontal outlet Γ out (see Figure 3). The aim is to nd the uid ow optimal domain Ω * which minimizes the dissipated energy by the forward Find Ω * solution to min where with |.| and V desired represents respectively the Lebesgue measure and the target volume. We recall that the performance of the Tesla valve is measured by diodicity Di which is known as the ratio of the pressure drop in backward direction to that in forward direction, which is equivalent to the ratio of dissipation of reverse and forward ows [15]: with (w f , p f ) and (wr , pr) are respectively the solution to the Navier Stokes system in the forward and the reverse ows. Then, diodicity can be maximized by minimizing forward dissipation while maximizing reverse dissipation. That is why our optimization problem is de ned with the diodicity Di as a constraint; Di > . Using the above de nitions, the optimization problem [15] for reconstructing Tesla valve can be expressed as * Objective: power dissipation of forward ow Navier Stokes equations for forward and backward directions.
We use the obtained theoretical results in 3.2.2 to solve (12).

. . The topology optimization process
To obtain the optimal domain, an iterative process is applied to construct a sequence of geometries (Ω k ) k≥ with Ω = Ω and Ω k+ = Ω k \ O k where O k is an obstacle inserted in Ω k . To de ne the obstacle location and size, we nd the function δj k de ned by (see Theorem 3.7) where − w k represents the velocity, solution to the Navier-Stokes problem in Ω k − u k is the adjoint state, solution to where Σ k = ∂(∪ k l= O l ) is the obstacle boundary inserted during the previous iterations. The optimization steps are summarized as: The Algorithm: 1. Initialization: Set Ω = Ω, and k = 2.
Repeat until |Ω k | ≤ V desired : (a) The topological sensitivity function: compute w k , solution to the non stationary Stokes problem (14) in Ω k , - compute v k , solution to the associated adjoint problem (15) in Ω k , -compute the term δJ k and deduce the function δj k (z), ∀z ∈ Ω k .
The obstacle to be inserted: The new domain: The stopping criteria is de ned by the natural optimality condition This algorithm is like a descent method where δj k represents the descent direction and |O k | = |Ω k \ Ω k+ | the step length. The parameter ρ * k is chosen to allow ρ −→ j(Ω k \ O k ρ ) to decrease as much as possible. The computation of ρ * k in (b) can be viewed as line search step. The numerical discretization of problems (14) and (15) is done by P − bubble/P nite element method [25]. The computation of the approximated solutions is achieved by the Uzawa's algorithm. The function δj k is computed piecewise constant over elements. Next, we will apply the proposed algorithm to reconstruct Tesla micro valve.

. . Reproducing the Tesla micro valve
We illustrate in this section the strengths of topology optimization method, namely the ability to nd optimal design using only information on boundary conditions and constraints without the need of initial design.
The considered design domain is the pentagon domain (see Figure 3). This problem example has already been studied by S. Lin and al. in [15] in the steady state regime using projection method.
For the forward direction, the inlet boundary velocity has a parabolic behavior (Re=100 relating to the inlet dimension). At the outlet boundary, the pressure is taken constant and no-slip condition is considered on the walls. For the backward ow direction, we reverse these boundary conditions. Besides, we prescribe solid regions close to the inlets/outlets to minimize the boundary e ect on the nal design solution. We illustrate the geometries obtained during the optimization process in Figure 4. The optimal domain is obtained after four iterations. It is nearly identical to literature [1,15] (see Figure 5). . .

Discussion
In the previous Figure, thanks to the topological gradient, we deduce an easy reconstruction of tesla valve. Now, we normalize the obtained tesla valve behavior by plotting the obtained forward and reverse ows respectively in gures 6(a) and 6(b). It is clear that the velocity eld is strongly di erent for the two cases. To study the obtained tesla valve performance, we calculate the diodicity. Using the energy view point expression of diodicity, the experimentally derived value is .
. In bibliography [26], the diodicity is well predicted using with N is the number of tesla valves and Re is the Reynolds number. Based on this expression, we found Di ∼ = .
which ensures an agreement between the obtained diodicity and the experimental one.

Mathematical analysis
This section deals with the proofs of Theorems 3.1, 3.3, 3.5 and 3.7.

. Proof of Theorem 3.1
Let Q be the pressure associated with the velocity W: where P j is the pressure associated with the velocity U j solution to (6). Setting the variation From (4) and (6), we can verify that (zε , pz ε ) is solution to The last boundary condition follows due to the fact that U j = −e j on ∂O. Moreover, since |w | L ( ,T;H (Ω)) < ν/k, then ε su ciently small, |w + W| L ( ,T;H (Ωz, ε )) ≤ α < ν/k.
Yet, we have • Estimate of boundary condition imposed on Γ: Let R > such that Oz, ε ⊂ B(z, R) and B(z, R) ⊂ Ω. Since z ∉ Ω R = Ω\B(z, R), the function x → E(x − z) belongs to C (Ω R ). By the trace theorem, we have Therefore, W L ( ,T;H / (Γ) is uniformly bounded with respect to ε.
• Estimate of boundary condition imposed on ∂Oz, ε: Using the theorem of trace and the smoothness of w in Oz, ε×] , T[, one can obtain Then, the rst boundary term on ∂Oz, ε satis es To estimate the last boundary term, we use that O contains the origin. Setting Or = O\B( , r) and Or, ε = z + εOr. Using the theorem of trace and the variable change x = z + εy, we obtain From the fact that y → E(y) is su ciently smooth in Or ⊂ R \ { }, the last quantity is uniformly bounded and then Finally, combining the above estimates, we obtain, ∃c > , independent of ε, such as zε L ( ,T;H (Ωz, ε )) ≤ −c log(ε) which ends the proof of Theorem 3.3.

. Asymptotic analysis
This section deals with the proofs of the Theorems presented in paragraphs 3.2 and 3.3. Using the assumption (A), where wε is extended by zero inside the domain Oz, ε.
Using Green formula and that wε = in Oε, it follows where u is the solution to the associated adjoint problem. From (4) Therefore, We begin by giving the estimate of the rst three terms in (33).
By the same arguments, we can estimate the two other terms in (33).
The shape function variation can be rewritten We are now ready to prove the established results in Theorems 3.5 and 3.7 and propositions 3.8 and 3.9.
The variation of the associated shape function j is given by ). Therefore

Conclusion
This paper deals with non-stationary Navier-Stokes topological optimization problem. In the theoretical part of this work, we have established a topological asymptotic formula describing the shape function variation related to a small Dirichlet geometric perturbation. The obtained theoretical results are exploited for building a topological optimization algorithm for solving the Tesla micro-valve optimization problem. We illustrate the strengths of this approach namely the ability to nd optimal design based only on boundary conditions and constraints information without the need of an initial design.