Shape minimization of the dissipated energy in dyadic trees

In this paper, we study the role of boundary conditions on the optimal shape of a dyadic tree in which flows a Newtonian fluid. Our optimization problem consists in finding the shape of the tree that minimizes the viscous energy dissipated by the fluid with a constrained volume, under the assumption that the total flow of the fluid is conserved throughout the structure. These hypotheses model situations where a fluid is transported from a source towards a 3D domain into which the transport network also spans. Such situations could be encountered in organs like for instance the lungs and the vascular networks. Two fluid regimes are studied: (i) low flow regime (Poiseuille) in trees with an arbitrary number of generations using a matricial approach and (ii) non linear flow regime (Navier-Stokes, moderate regime with a Reynolds number 100) in trees of two generations using shape derivatives in an augmented Lagrangian algorithm coupled with a 2D/3D finite elements code to solve Navier-Stokes equations. It relies on the study of a finite dimensional optimization problem in the case (i) and on a standard shape optimization problem in the case (ii). We show that the behaviours of both regimes are very similar and that the optimal shape is highly dependent on the boundary conditions of the fluid applied at the leaves of the tree.


Introduction and motivations
Tree structures are very common means to transport a product between two regions of different scales. A fluid acting as a transporter often flows in such structures. The circulation of the fluid in such geometries can dissipate a lot of energy by viscous effects and the question of the optimization of the tree geometry arises. Important applications of this problem exist, like in industry, in the study of river basins, in water treatment, in medicine. Such structures are also often encountered in Biology and are the result of evolution. One can see roughly natural selection as an optimization process which adapts the organisms to their environment. Thus, a major question arises for biological systems: what are they optimized for? We know that an effect of natural selection is to minimize some complex cost functions relatively to a wide range of parameters. Although most of the time unknown, these parameters have probably various influences in term of amplitude, some stronger than others. Thus, if we are able, typically through a modeling work, to give hypotheses on what are the most influential parameters, to isolate them in a model and to determine their role if they were alone, then a simple comparison between the results of the model and the real biological system could indicate whether their role was truly important or not. In the case of biological networks such as lungs or vascular network, two cost parameters arise naturally: the viscous dissipated energy of the fluid in the network and the volume of the network. Indeed, these organs have to deal with energy dissipation due to air or blood circulation [34] and since they span in the geometry they have to feed, they cannot use too much volume. These organs can be subject to dysfunctions which are often consequences of an increase of their hydrodynamical resistance and thus of a loss of efficiency of the geometrical structures as transport systems. Hence, dysfunctions like asthma or heart attacks are often linked to an increase of the hydrodynamic resistance (energy cost spent for the circulation of the fluid) of the structure. Thus, researching the optimal shapes of tree structures could bring useful information. Previous studies have been made on this topic in the past like in [16,2,33,20,3] each on particular situations and applications.
In this frame, the goal of this work is to determine the shapes of dyadic trees that would minimize the viscous energy of a Newtonian fluid under a volume constraint on the tree. As written upwards, such structures are good candidates for the modeling of mammals bronchial trees [20] and we will focus most particularly on this application. We study this problem for two regimes of flow. We begin with low flow regime (Poiseuille flow) using a matricial formulation of the problem as in [22,21]. We also study the non linear Navier-Stokes flow at a moderate Reynolds number (∼ 100) which is however sufficient to exhibit inertial effects around the bifurcation. In this second case, we use a numerical shape minimization method based on shape derivatives [13]. We show that for both regime the optimal structure depends on the fluid boundary conditions that are applied at tree root and leaves. Indeed, under the assumption that the total viscous flow is constant throughout the tree, the optimal shape is very different according to the boundary conditions imposed at the leaves: Dirichlet conditions or strictly identical Neumann conditions lead to an optimal tree with particular relationships between the flow in a branch and its diameter; non identical Neumann conditions lead to a degenerated tree reduced to a tube and only one leaf remains accessible to the fluid from the root. Moreover the numerical simulations in the non linear case give precisely the geometry of the bifurcation. We focus here on the 3D case, however it is easy to see, with very few changes in the reasoning, that our results would also hold for the 2D case.

Terminology and notations
In the following, we will call inlet of a dyadic tree either the open surface of the root of the tree or the root of the tree itself, depending on the context. Similarly, we will call outlets of a tree either the open surfaces of its leaves or its leaves themselves, also depending on the context. Mainly, we will refer to the open surfaces if we speak of boundary conditions and to the branches in the other cases. Note that this terminology does not mean necessarily that fluid is going in the tree through the inlet and out of the tree through the outlet.
We will use the following notations throughout this section: [x], with x ∈ R the integer part of x B ⊤ , where B is a matrix the transpose of B (e 1 , . . . , e n ), with n ∈ N * the canonical basis of R x, y , where x and y are two vectors with same length the euclidean inner product .
the euclidian norm, induced by the inner product ., . diagu, where u = (u i ) i∈ 1,n ∈ R n×1 the diagonal matrix D = (d i,j ) 1≤i,j≤n such that d i,i = u i for all i ∈ 1, n .
In the following, a family of real numbers (ν ε ) ε∈R * + is assimilated to a sequence since for each n ∈ N * , one can choose ε = 1 n .
Moreover, let us define the direct sum of two matrices.
Definition 1. Let m and n be two nonzero integers. Let A 1 ∈ R n×n and A 2 ∈ R m×m be two matrices. The direct sum of A 1 and A 2 is the matrix M in R n+m,n+m defined by blocks by

Models
We introduce here all the models used in this article. For the sake of clarity, all proofs of this section have been regrouped in Appendix A.

Poiseuille's law
We consider a viscous incompressible fluid whose dynamic viscosity is µ > 0. It flows in a steady and laminar state through a cylindrical rigid pipe whose length is L > 0 and radius is R > 0. We impose a no-slip condition on the lateral boundary which means that the fluid "sticks" to the wall. We refer the inlet to 0 and the outlet to 1. Pressures (p 0 > 0 and p 1 > 0) at its openings are supposed to be uniform all over the section. The volumetric flow rate Φ is chosen positive if the flow goes from section 0 to section 1. The fluid behavior is ruled by the Navier-Stokes equations in the cylinder. The solution of these equations with the previous conditions is characterized with velocity profiles that are parabolic on each section and a pressure that is constant on each section and decrease linearly along the axis of the pipe. This type of flow is known to verify Poiseuille's law. This law boils down to a linear relationship between the pressure drop p 0 − p 1 and the volumetric flow rate Φ through the pipe, that is Therefore, if the pressure drop p 0 − p 1 is positive, the flow goes from section 0 to section 1.
There is an exact correspondence between Poiseuille's law and Ohm's law if we match pressure drop with potential difference and volumetric flow rate with electric current. Thus by analogy, the proportionality constant r is called a hydrodynamic resistance. Finally, since the flow is incompressible, the flow rate is conserved through the pipe.

Dyadic trees
In this section, we will define the different mathematical tools needed in the sequel to manipulate tree structures in which flows a Poiseuille's fluid. In particular, we give a matricial relationship between the flow at the inlet of a dyadic tree and the pressures at the outlets.
We consider the flow of an incompressible and viscous fluid whose dynamic viscosity is µ through a finite dyadic tree of N + 1 generations (N ∈ N, N ≥ 1). We recall that a new generation of this tree occurs when a bifurcation is created. Hence the root branch corresponds to generation 1 and the branches at the leaves of the tree corresponds to generation N + 1. Consequently, the tree has 2 N outlets and 2 N +1 − 1 branches.
We also define the level as a number associated to a generation, equal to 1 for the second generation and increased by 1 at each generation. Therefore, a tree with N + 1 generations has N levels.
We do the distinction between the generation and the level in the tree, because it makes the indexation of the different variables of the tree easier. The levels will be denoted in all this section by the index letter i.
We assume that our tree is composed of connected rigid cylindrical pipes in which the fluid obeys to Poiseuille's law. We call Φ > 0 the volumetric flow rate that enters the tree. Since the flow is incompressible, the flow rate is conserved through the pipes and at bifurcations.
The sets B N,i of couples of indexes that locate each branch of a given level i (i ∈ 1, N ) is Hence, i represents the level and j the position of the branch at this level. Thus, the set B N of all indexes locating the pipes for the overall tree is This set does not include the root branch of the tree. This branch has a radius R 0 > 0, a length L 0 > 0, the pressure at its inlet is p 0 > 0 (pressure at the inlet of the tree) and is p 1 > 0 at its outlet. Thus, its hydrodynamic resistance is r 0 = 8µL 0 πR 4 0 . According to Poiseuille's law, one has For a pipe whose location in the tree is given by the couple (i, j) ∈ B N , we denote R i,j > 0 its radius and L i,j > 0 its length. Therefore, the hydrodynamic resistance of this pipe is We use q i,j and p i,j > 0 to define respectively, the volumetric flow rate through this pipe and the pressure at its outlet. The flow rate q i,j is chosen positive if the fluid in a pipe flows towards the pipes with a higher generation index.
We consider now a pipe whose index is (i, j) ∈ i∈ 1,N −1 B N,i . Since the tree is dyadic, the hydrodynamic resistances of its (two) daughter pipes are r i+1,2j−1 and r i+1,2j , which belong to the (i + 1)-th level. Then, we define the reduction ratios x i+1,2j and x i+1,2j+1 , that represents the change in the geometry of the pipes between the levels i and i + 1, by (3) Moreover, we assume an identical reduction ratio for the radius and the length between a mother branch and its daughter, thus (4) Figure 1 shows an example of a dyadic tree and of our notations in the special case where N is equal to 2. Our goal is to establish a relationship between the 2 N pressures and the 2 N volumetric flow rates at the outlets of the tree. To go further, we need to be able to follow the fluid through paths in the tree. Therefore,we define the notions of path and subpath in the tree.
Definition 2. Notion of path and subpath.
is the set of couples of indexes of the i branches needed to link the root branch denoted by 0 and the branch located by (i, j) in the tree. It includes the branch referred to by (i, j) but not the root branch. More precisely, where (m k ) 1≤k≤i denotes the sequence of positive integers defined by 2. Let (i, j) ∈ B N and s ∈ 0, i . For a given path Π 0→(i,j) , the subpath Π 0→(i,j) (s) is the set of couples of indexes of the m branches needed to link the root branch and a branch located at the s-th level following a part of the path Π 0→(i,j) . More precisely, the subpath Π 0→(i,j) (s) is the subset of Π 0→(i,j) defined by We use this definition to do a change of variable. It will be very useful in the second part of this section devoted to the optimization of a given criterion with respect to the geometry of the tree represented by the 2 N +1 − 2 variables {x i,j } (i,j)∈B N . From now on, we replace the variables x i,j by the new ones ξ i,j defined by We will impose another constraint on the geometry: lengths and radii of the tree are assumed to decrease as we go along its levels. More precisely, It has to be noticed that the map {x i,j } (i,j)∈B N −→ {ξ i,j } (i,j)∈B N is obviously a C 1 -diffeomorphism on the set of strictly positive real numbers so that it defines a change of variable. Moreover, r 0 ξ i,j represents the hydrodynamic resistance r i,j of the pipe denoted by (i, j). For instance, in the case N = 2 (see figure 1), the path linking the inlet to the first outlet is Π 0→(2,1) = {(1, 1); (2, 1)} and the geometric variables of theses branches are , ξ 1,1 = x 1,1 , and x 2,1 = r 0 Since the tree is dyadic, it is usual to compute the total volume V of the tree by summing the volume of each cylindrical branch, i.e.
Notice that this volume is an approximation of that of the real tree, since it takes into account only the volumes of the cylinders composing the tree and not the volumes of the bifurcations.

Now, let us define
• the vector p containing all the pressures at the outlet of the tree, i.e.
• the vector q containing all the volumetric flow rates at the leaves of the tree, i.e.
Definition 3. Given two positive integers i and j and their binary expressions Thanks to Poiseuille's law, we are able to establish a linear relationship between the vectors p and q.
where • the symmetric matrix A N (ξ) ∈ R 2 N ×2 N is called resistance matrix of the tree and is defined by • u N ∈ R 2 N ×1 denotes the ones vector i.e. u N = (1, . . . , 1) ⊤ .
is the subpath corresponding to the intersection of the two paths Π 0→(N,i) and Π 0→(N,j) .
This proposition is proved in appendix A.1. A close result has already been proved in [22].
The quantity E(q, ξ) = q ⊤ A N (ξ)q (12) corresponds to the total viscous energy dissipated by the fluid in the tree during one second, see [22,21], i.e.
Using the fact that the flow rate is conserved all along the pipe, we rewrite the dissipated energy as As a consequence, A N (ξ) is a symmetric positive definite matrix. It justifies the following proposition.
Example 1. So as to have an idea of the structure of the matrix A N (ξ), let us use an example of a tree with N = 2 levels (thus, with N + 1 = 3 generations, 2 N −1 = 4 outlets and 2 N − 1 = 7 branches). Its resistance matrix A 2 (ξ) ∈ R 4×4 is defined as follows: The tree corresponding to this example is drawn on the figure 1.

Boundary conditions
For biological systems like the vascular network or the lung, it is very difficult to know what are the correct boundary conditions since there exists complex feedback loops that are able to modify the forces applied by muscles relatively to the behaviour of physiological values. Typically, lungs ventilation is controlled not only by carbon dioxide and oxygen concentration in blood but also by blood pH and by mechanical sensors distributed along the tree [28]. Concerning boundary conditions in the frame of fluid mechanics, one can refer to [5] and [17] and in the context of blood flow modeling to [27]. We will investigate the two major types of boundary conditions at the openings of the tree. We will either impose a flow (or parabolic velocity profile) which corresponds to a quantity of fluid that would circulate in the branch, or impose a pressure which corresponds to a force that would be applied to the opening of the branch. If flow rates are imposed at the outlets of the tree, then the pressure will be imposed at the inlet and vice versa so that the problem is well posed. Indeed, the two cases studied in this work are • 1st case: pressure is imposed at the inlet and flow rates at the outlets.
• 2nd case: flow rate is imposed at the inlet and pressures at the outlets.
In the 2nd case, the flows at outlets are not directly accessible. Thus, the following proposition gives the existence of q as the solution of a linear system. Proposition 3. Assume that the pressures at the outlets of the tree are fixed and the flow in the tree root equal to Φ. Then, the vector q ∈ R 2 N ×1 is the unique solution of the linear system where and

The optimization problems
In this section, we will study the two finite-dimensional constrained optimization problems corresponding to two types of boundary conditions: • 1st case: pressure is imposed at the inlet and flows at the outlets, • 2nd case: flow is imposed at the inlet and pressures at the outlets.
We consider the same rigid dyadic tree as in section 2 with N + 1 generations (N ∈ N, N ≥ 1) through which flows the fluid introduced previously. The tree is characterized by its geometry ξ and the volumetric flow rates q (q ∈ R 2 N ×1 ) at its outlet. We recall that, at outlets, pressures are linked to flow rates with the relation p = A N (ξ)q.
We want to minimize the total viscous dissipated energy E defined by (12) with respect to the pair (q, ξ) under the constraints • that q is the flow vector at outlets of the tree, when the Poiseuille's model is considered. This condition can always be treated as a constraint, it takes however different forms depending on the case considered.
• and that the total volume V of the tree is constrained, which from (9), writes To make the admissible set of pairs (q, ξ) clear, we define the intermediate set A Λ as (8) and (16) .
3.1 1st case: pressures at the inlet and flow rates at the outlets are known The first case is much simpler than the second one (presented in section 3.2). The pressure is imposed at the inlet along with the flows vector q at the outlets, thus the optimization problem is where q ∈ R 2 N ×1 + is assumed fixed. In this case, since the volumetric flow rate is conserved along the tree, every value of the intermediate flow rate q i,j is known from the values of q N,j , j ∈ 1, 2 N .
The existence of solutions for the problem (P 1 ) is classical. Indeed, if one of the optimization variables goes to zero then the energy tends to a positive infinite value. It is thus possible to come down to minimize E(q, ·) on a compact set, which yields the existence of a solution. Moreover, A Λ is a convex set and E(q, ·) is strictly convex, this ensures the uniqueness of the minimizer for the problem (P 1 ).
We will now write the first order optimality conditions: we denote by ξ * the minimizer of the problem (P 1 ). There exists a Lagrange multiplier λ ∈ R such that Using the volume constraint, we immediately obtain the following expression for the minimizer ξ * : The following proposition summarizes the conclusions for the first case.
Proposition 4. The problem (P 1 ) has a unique solution ξ * given by

2nd
Case: flow rate at the inlet and the pressures at the outlets are known The volumetric flow rate Φ > 0 which enters the overall tree through the root branch and the pressures at the outlets of the tree (i.e. the vector p, p ∈ (R * + ) 2 N ×1 ) are assumed to be fixed. According to Proposition 3, the constraint q = (M N (ξ)) −1 b N comes down to impose an incompressible and Poiseuille-like flow through the tree. Now, let us define the set U Λ of admissible pairs (q, ξ) as The resultant optimization problem writes We claim that generically with respect to the vector p, the problem (P 2 ) has no solution. Nevertheless, in the degenerate situation where all the pressures at the outlet of the tree are equal, the problem (P 2 ) has a unique solution.
We will start with the case where (P 2 ) has a solution. The following theorem gives a characterization of this situation. Theorem 1. The problem (P 2 ) has a solution if, and only if there exists π * ∈ R such that This theorem emphasizes the fact that if two pressures at the outlets are different, then the problem (P 2 ) has no solution. In this case, we are able to exhibit a minimizing sequence. It has to be noticed that the value of the infimum does not differ according to the case considered.
Theorem 2. Assume that at least two pressures differ at the outlets of the tree, i.e.
The limit case "q 0 = (Φ, 0, . . . , 0) ⊤ " corresponds to the case where the fluid exits the tree only by the outlet denoted by (N, 1). The use of some symmetry properties in the structure of the objective function E yields that we would easily get another minimizer by choosing any other outlet as main exit of the fluid, which would correspond to the closure of every path linking the root branch to the outlet except Π 0→(N,i) , for any i ∈ 1, 2 N . Note that the symmetry property does not hold for the pressures at outlets. However, when we take the limit in the minimizing sequence, since the total flow in the tree is imposed, the pressure drop between the root and the outlet that remains open is the same whatever the outlet chosen. Thus, the pressure at root will depend on the pressure imposed at the outlet. Hence, introducing the matrix M N (ξ) that does not depend on the root pressure is well adapted.
The optimal shape for the case with each pressures equal at outlets is obviously unstable relatively to the pressures: a slight change in one of them will shift the optimal shape from a tree-like structure to a pipe-like structure. From the modeling point of view, it means that a tree-like structure verifying the hypothesis of this section should not be encountered in nature as a result of an evolution process such as natural selection. Indeed a perfect adjustment of pressure at outlets is very difficult to assure and maintain.
The proofs of the theorems will be decomposed into two steps, whose details can be found in the next section: 1. determination of a lower bound for the total viscous dissipated energy, 2. construction of the minimizing sequence (q ε , ξ ε ) ε>0 .
Notice that, in the rest of the paper and besides Section 4.1, the notation q points out the vector composed only of the flow rates at the outlets of the tree. Let r 0 > 0, Φ > 0 and Λ > 1. Let (Q) be the finite-dimensional constrained optimization problem where E(q, ξ) is defined by and the sets of constraints are Moreover, the value of the minimum is equal to Note that there is no reason to have uniqueness of the solution of Problem (Q).
Remark 2. The expression of the necessarily first-order optimality conditions in the proof of Proposition 5 leads easily to conclude that Problem (Q) has an infinite number of minimizers.
Moreover, an interesting minimizing sequence (q ε , ξ ε ) ε>0 in C 1 × C 2 , of E for our problem is given by Physically, it corresponds to the case in which the fluid exits the tree only by the outlet denoted by (N, 1).

Proof of Theorem 1
As a preliminary to the proof, let us deal with the qualification issue for Problem (P 2 ). In fact, any elements of the set U Λ verifies the constraint qualifications. Indeed, thanks to Proposition 3, one can write q = M N (ξ) −1 b N and furthermore, it is obvious that any element of the set A Λ verifies the constraint qualifications. Combining the two previous remarks, the set U Λ inherits from this property . Let us assume that Problem (P 2 ) has a solution (q * , ξ * ). In the proof of Theorem 2 presented in Section 4.3, it is proved without any additional assumption than those of Theorem 1, that the sequence (q ε , ξ ε ) ε>0 defined by (22) is a minimizing sequence for Problem (P 2 ), whatever values of pressures at the outlets. Hence it follows that Now, recall that, by Proposition 5, any solution of Problem (Q) realizes the minimum r 0 Φ 2 1 + N 2 Λ−1 and verifies the necessary optimality conditions (25). Since E(q * , ξ * ) = r 0 Φ 2 1 + N 2 Λ−1 and (q * , ξ * ) ∈ U Λ ⊂ C 1 × C 2 , the pair (q * , ξ * ) belongs to the set of minimizers of Problem (Q) and verifies therefore, the necessary optimality conditions (25), namely Let us denote as previously by p 0 the pressure at the inlet of the tree. According to Poiseuille's law and reasoning by induction, we obtain In particular, it implies that Conversely, let us prove that, if there exists π * ∈ R such that p = π * u N , then, Problem (P 2 ) has a solution. In this case, we are able to exhibit an element (q * , ξ * ) which belongs to the admissible set U Λ . For instance, let us choose Then, it is easy to make the calculation of the flow throughout the tree. Indeed Iterating this reasoning proves that q * is proportional to u N and finally, by induction, we obtain A direct calculation shows hence that

Proof of Theorem 2
The fact that Problem (P 2 ) has no solution when assumption (21) holds is a direct consequence of Theorem 1. In this section, we prove that the sequence exhibited in (22) is a minimizing sequence for Problem (P 2 ), in other words, where q(ξ ε ) is the unique solution of the linear system M N (ξ ε )q = b N (see Proposition 3). Indeed, since U Λ ⊂ C 1 × C 2 and since the sequence (q(ξ ε ), ξ ε ) ε>0 has been constructed so that its elements belong to U Λ , one has for ε sufficiently small, Hence, proving (33) will yield at the same time that (q(ξ ε ), ξ ε ) ε>0 is a minimizing sequence for Problem (P 2 ) and that the infimum is equal to r 0 Φ 2 1 + N 2 Λ−1 .
Let us denote for the sake of clarity, q(ξ ε ) by q ε , A N (ξ ε ) by A N (ε) and M N (ξ ε ) by M N (ε). Now, recall that, according to Proposition 3, the vector q ε is completely characterized by the system Let ε > 0. The matrix A N (ε) admits the decomposition where U N , A 1 N and A 2 N are three matrices independent of ε with same size as A N (ε), such that • The first part, U N is the matrix of size 2 N × 2 N whose coefficients are all equal to 1. This term comes from the root branch that add a resistance r 0 to the resistance of any pathway in the tree.
• The second part, r 0 ε A 1 N , corresponds to pathways consisting only in branches (i, j) such that ξ (i,j) = ε (i.e. whose diameter will go to 0 with ε), namely all pathways in the tree except those included in the pathway going from the root to the outlet (N, 1). A 1 N is more precisely characterized in Lemma 1 below.
• The third part N corresponds to the pathways that stay open when ε goes to 0.
Let us define A N (ε) = εA N (ε). Then, we have Therefore, the vector q ε is completely characterized by the system Thus it is necessary to study more closely the matrix A 1 N that will give the behaviour of q ε when ε goes to 0.
is a resistance matrix of a tree of p + 1 generations whose branches have all the same resistance equal to 1.
Proof. From its definition, A 1 N corresponds to the resistance matrix of a tree with N +1 generations where the resistances of the branches on the path from root to outlet (N, 1) are all 0 and the resistance of the other branches are all 1. The numbers next to the branches represent their resistance. The matrix A 1 N is the direct sum of the resistance matrices of the subtrees built up by branches with non zero resistance.

Now, system (37) is equivalent to
Proof. Let us show that the family ( We consider some real numbers (λ i ) i=1...2 N −1 and µ such that Let us recall that, from Lemma 1, B p where B p is a resistance matrix and is thus invertible. Consequently, the range of A 1 N is Span (e 2 , . . . , e 2 N ), which implies that, in equation (39) the component along e 1 is reduced to µ = 0. Now, the restriction of A 1 N on its range is a one-to-one correspondence as soon as the projection of the family (v i ) i=1...2 N −1 on the range of A 1 N is a basis of the range of A 1 N , which is clearly true in our case. Then, the family ( Then, since we have b N = εγ + Φe 2 N , with γ ∈ Span e 1 , . . . , Since lim ε→0 (εA N (ε)) = A 1 N and lim ε→0 M N (ε) = M N , then one successively has which concludes the proof.

Case of a fluid driven by Navier-Stokes equations: some numerical results
As stated in the introduction, we study in this section an infinite dimensional optimization problem. Now the unknown is the whole shape of the dyadic tree and we do not anymore neglect the connections at bifurcations nor do we constrain the branches to remain cylindrical. Moreover, we will now use the full non linear Navier-Stokes equations and still minimize the dissipated energy of the fluid. Our goal is to compare the shapes obtained in this case by numerical means to that obtained theoretically for Poiseuille's laws in the previous sections. Nevertheless, we will not be able to catch the behaviour of the case with identical pressures imposed at each outlets since it is unstable. It is not possible in our numerical simulations to make the pressure exactly equal because of the rounding errors and the local meshes. Only a dedicated algorithm that would include the symmetry and prevent the closing of branches would be able to catch this behaviour. Consequently, from the theoretical study, we expect to observe the closing of every branches except one. We focus here on the 3D case but the 2D case can be easily adapted from this 3D study and numerical results will be presented both in 2D and 3D.
The partial differential equations describing the behaviour of the fluid and the boundary conditions is now more general. In particular, this model is a convenient choice for the modeling of the bronchial tree (see for instance [18,19,20,22]).
Let us recall some classical definitions in incompressible Fluid Mechanics.
Definition 4. Let u be a smooth vector field of R 3 (for instance C 1 ). We define 1. the stretching tensor of u (symmetric part of the gradient tensor): 2. the doubly contracted product of two stretching tensors ε(u) and ε(v): 3. the stress tensor of (u, p), where u is a smooth vector field of R 3 representing the velocity of a fluid whose viscosity is µ > 0 and p a function representing the pressure defined on R 3 : where I 3 is the identity tensor of R 3 .

The shape optimization problem
We now give some precisions on the frame of our study. Let Ω be a generic three dimensional Y -shaped domain. We will denote by ∂Ω the boundary of Ω. In the sequel, we will assume that the inlet E of Ω is a disk and that the outlet S consists in two identical disks. We decompose the boundary of Ω as the disjoint union ∂Ω = E ∪ Γ ∪ S. Γ, the lateral boundary, is the main unknown or the shape we want to determine. In the following, we will thus consider that E and S are fixed part of the boundary. We assume that a fluid driven by the Navier-Stokes equations flows in Ω. In particular, this model is convenient to represent the upper part of the bronchial tree since the velocity of the air at the beginning of the trachea is important enough to consider that the flow is inertial or even turbulent (the Reynolds number in trachea is about 1000 at rest). Then, the velocity u : Ω → R 3 and the pressure p : Ω → R, are solutions of the Navier-Stokes system where • n denotes the outward-pointing unit normal vector at a given point of the boundary ∂Ω, • u 0 is a parabolic velocity profile (i.e. a Poiseuille's flow is imposed at the inlet E), that is where c is a negative constant so that the flow is ingoing, and R the radius of the inlet.
Remark 3. Let us clarify the choice of the boundary conditions for this model.
• The condition u = 0 on Γ is the so-called no-slip boundary condition and means that the fluid "sticks" to the wall.
• Practically speaking, we will impose h = −p 0 n at the outlet S. This condition comes more or less from the assumption that the pressure is the sole force acting on S and that it is known and equal to p 0 . Drawing a parallel with the modeling of the bronchial tree, these conditions simulate muscles applying the same force per unit of surface for each outlet (−p 0 n) and thus using the same energy for each (see [22]). Similar boundary conditions are more detailed in [5, chapter 5] and in [6,7].
Notice that the classical theory of Navier-Stokes equations (see [11,32]) gives information on the existence and uniqueness of the solution of System (42): Let Ω be a bounded Lipschitz domain of R 3 . Let us assume that u 0 belongs to the Sobolev space (H 3/2 (E)) 3 and h belongs to (H 1/2 (S)) 3 . There exists A > 0 such that if the viscosity µ is larger than A, then Problem (42) has a unique solution (u, p) ∈ (H 1 (Ω)) 3 × L 2 (Ω).
Let us now introduce the shape optimization problem we want to solve. The objective functional J(Ω) is the energy dissipated by the fluid, i.e.
To make the statement of the optimization problem precise, we need to define a class of admissible shapes. As classically in shape optimization, we fix the measure V 0 > 0 of Ω. Let us introduce O V 0 = Ω bounded and simply connected domain in R 3 containing E and S, such that meas(Ω) = V 0 .
(44) The shape optimization problem writes The question of knowing if Problem (45) has or not a solution is still an open problem. Nevertheless, it is possible to show that a problem, very close to Problem (45) but a little bit constrained, has a solution.
Restricting the set of admissible shapes is a very common approach in shape optimization, since these problems are often ill-posed (see for instance [1,13]). A very close existence theorem is announced in [14] and proved in [15], considering instead of the set O V 0 , a set of domains verifying an ε-cone property, which yields some kind of uniform regularity (see [8,13,9] for some reminders on the ε-cone property and its consequences on the existence of optimal shapes).

Computation of the shape derivative of J
The classical algorithm we will recall in Appendix B is based on the use of the well known shape derivative. We recall here the expression of such a derivative. The details of its calculation are given in [15].
It is possible to define an adjoint state for System (42) that is useful to write the shape derivative in a very usual form (see [13,Theorem 5.9.2]). Let (v, q) be solution (in the case where it exists) of the linearized Navier-Stokes system The following existence and uniqueness result is established in [15, Proposition 3.1].

Proposition 6.
Let Ω be a bounded Lipschitz domain of R 3 . There exists A > 0 such that, if the viscosity µ is larger than A, then the problem (46) has a unique solution (v, q). Moreover, this solution belongs at least to (H 1 (Ω)) 3 × L 2 (Ω).
Notice that the restriction in the size of the viscosity µ ensures the existence, the uniqueness and a sufficient regularity of u, solution of the Navier-Stokes equation 42 (see theorem 3).
We are now able to define the derivative of J with respect to the domain. Let us consider a regular vector field V : R 3 → R 3 with compact support which does not meet neither E nor S. For small t, we define Ω t = (I + tV )Ω, the image of Ω by a perturbation of identity and f (t) = J(Ω t ). We recall that the shape derivative of J at Ω with respect to V is f ′ (0). We will denote it by dJ(Ω; V ). To compute it, we first need to compute the derivative of the state equation. We use here the classical results on shape derivatives as in [13], [24], [31]. The derivative of (u, p) is the solution of the linear system where u is, under assumption that µ is large enough, the unique solution of (42). Similar formula in the context of Shape Optimization in Fluid Mechanics are established for instance in [23,25,29,30]. Now, we have (see [13], [31]) Using the adjoint state (46), it is possible to rewrite dJ(Ω, V ) in a more workable form (see [15,Proposition 3.2]): We will now introduce the shape gradient of the criterion J as a functional defined on the boundary such that dJ(Ω; V ) = Γ ∇J(Ω) · V ds, in other words Thanks to this expression of the shape gradient, we are now able to implement an augmented Lagrangian algorithm in order to solve Problem (45). The algorithm is classical and is described in Appendix B.

Numerical results
In this section, we present 2D and 3D numerical simulations using the augmented Lagrangian algorithm described in Appendix B based on the shape gradient (50). They have been implemented with the script of the software Comsol Multiphysics. The 2D algorithm is an easy adaptation of the 3D algorithm.
The Navier-Stokes systems and the adjoint state are solved with a direct finite elements method using Lagrange elements P 1 for pressures and P 2 for velocities. The displacements of the mesh are P 1. The multifrontal package (UMFPACK) is used to solve the resultant linear system and a modified Newton like method is used to treat the non linear term. At each iteration k, each node of the geometry is perturbed by the discretized operator (I + ε k d k )(Ω k ) (classical Arbitrary Lagrangian Eulerian method) and it is necessary to remesh the geometry when the displacement of the mesh becomes too large.
The Reynolds number in the following computations is 100 which is not only large enough to observe inertial effects near the bifurcation, but small enough to let the computations run in a reasonable amount of time. The effect of inertia on flow distribution in the bifurcation can be seen on Figure 3 where the image on the left represents a non inertial flow (Reynolds 0) and the image on the right an inertial flow with Reynolds 100. We normalized the applications J(·) and V (·) in such a way that J(Ω 0 ) = 1 and V (Ω 0 ) = 1.

1st Case: Neumann conditions are imposed at the inlet and Dirichlet conditions are imposed at the outlets
In this section, the initial geometry Ω 0 is a bifurcation whose branches are identical except for the length of the mother branch which is 1.5 the length of the daughter branches, see Figure 4 (2D, left image) and Figure 6 (3D, left image). The length of the mother branch has been chosen longer to avoid that the flow near the bifurcation interacts too much with the inlet flow since the algorithm tends to shorten the mother branch. Practically, a mother branch too short leads to convergence problem. The angle between two nearby branches is π/3. The boundary conditions correspond to that of Case 3.1 adapted to the non-linear regime. Thus a parabolic velocity profile is imposed at each outlets (Dirichlet condition) and Neumann boundary conditions σ(p, u)n = −p 0 n are imposed at the inlet, which is "formally" interpreted as a pressure constraint since viscous effects are negligible. The pressure imposed at the inlet is 0. Figures 5 (2D) and 7 (3D) show the convergence of the different quantities of the problem: the Lagrange multiplier (upper left), the Lagrangian function (upper right), the viscous energy (lower left) and the volume (lower right). Except the Lagrange multiplier, these quantities have been normalized. As expected, the curves are oscillating around their asymptotic values since verifying the volume constraint and minimizing the criterion are advantaged one after the other by the algorithm. This "advantage" is controlled by the Lagrange multiplier. Hence when it is too large the volume constraint is stronger and the algorithm authorized an increase of the criterion; if it is too small then the algorithm let loose the constraint and decrease the criterion. The amplitude of these oscillations decreases with the iterations.
In 2D, the mesh consists in 9296 triangle elements and no remeshing was needed to ensure convergence. Convergence was reached for 1000 iterations of the algorithm in about 6 hours using 2 cores of a Xeon processor (2.33 GHz). The viscous energy dissipated in the bifurcation, which is the quantity minimized by our algorithm, has been reduced by 20.4% in the final geometry in comparison with its value in the initial geometry. Volume precision at convergence was smaller than 2 × 10 −6 . In 3D, the mesh consists in 17959 tetrahedral elements and no remeshing was needed either. Convergence was reached for 1700 steps. Computation time was about 46 hours using 2 cores of a Xeon processor (2.33 GHz). The viscous energy dissipated in the bifurcation has been reduced by 22.5% between the initial and final geometry. The volume constraint at the end of the algorithm was smaller than 3 × 10 −5 . In this section, the initial geometry Ω 0 is a bifurcation whose branches are all identical, see Figure  8 (2D, left image) and 11 (3D, left image). The angle between two nearby branches is π/3. The boundary conditions correspond to that of Case 3.2 adapted to the non-linear regime. Thus a parabolic velocity profile is imposed at the inlet (Dirichlet condition) and Neumann boundary conditions σ(p, u)n = −p 0 n are imposed at the outlets. The pressures imposed at the two outlets are slightly different: 0 P a on one outlet and 2 × 10 −4 P a on the other.
As for the first case, the curves plotted on figure 10 (2D) and 12 (3D) represent the convergence of the different quantities of the problem: the Lagrange multiplier (upper left), the Lagrangian function (upper right), the viscous energy (lower left) and the volume (lower right). Except the  Lagrange multiplier, these quantities are still normalized. As before, they are oscillating around their convergence value (see 1st case).
In 2D, the mesh consists in 8624 triangle elements and no remeshing was needed. The convergence was achieved in 2800 steps for a computing time of 8 hours on 2 cores of a Xeon processor (2.33 GHz). In the final geometry the viscous dissipated energy is reduced by 54.5% relatively to the initial bifurcation. The precision of the volume constraint is smaller than 7 × 10 −5 .
In 3D, convergence was achieved in 7000 steps. A remeshing was performed at the step 4500 (represented by the dotted vertical line on figure 11). The initial mesh consisted in 10044 tetrahedral elements and the second mesh was finer with 37167 tetrahedral elements. The computation was much slower after the remeshing because of the increased number of elements. The total time needed by the algorithm to converge was of 259 hours on 2 cores of a Xeon processor (2.33 GHz). The viscous energy dissipated in the final geometry is reduced by 53.6% relatively to the viscous energy dissipated in the initial geometry. The volume constraint precision at convergence was smaller than 7 × 10 −5 .
Similarly than the result found theoretically at low regime in Section 3.2, we observe the closing of one branch. This indicates that this result should be more generally true and that it can probably be extended to non linear regime. A zoom of the 2D bifurcation is plotted on Figure 9, the arrows represents the normalized velocities and a region of fluid recirculation appears at the beginning of the closing branch. The colors on Figure 9 represents the pressure which is almost constant along the branch that is closing (up). This induces a very small flow inside it (the flow and the branch diameter are decreasing together when the optimization algorithm is progressing).

Conclusion
In this work, we show that fluid boundary conditions play an important role for the determination of the optimal shape of a tree in term of viscous dissipation. Indeed, the optimal shape associated to pressure conditions at leaves is a simple pipe while the optimal shape associated to flow conditions  at leaves is a tree. Moreover we have shown that these results hold for both Poiseuille's regime and a regime with a Reynolds number 100, which, although moderate, is large enough to exhibit inertial effects near the bifurcation.
Boundary conditions can be seen as constraints in the system and a pressure constraint is very different of a flow constraint. Let us consider a pipe in which flows a fluid in Poiseuille's regime. If we call R the hydrodynamic resistance of the pipe, the dissipated energy is J = RΦ 2 = (∆p) 2 /R (Φ is the flow, ∆p is the pressure drop). Then minimizing J relatively to the geometry of the pipe depends on the constraint: Situation (i): if the flow Φ is given then minimizing J is equivalent to minimizing the resistance R and the minimizer corresponds to ∆p = 0 (pipe radius goes to infinity).
Situation (ii): if the pressure drop ∆p is given then minimizing J is equivalent to maximizing the  resistance R and the minimizer corresponds to Φ = 0 (no flow in the pipe, the pipe radius goes to 0).
Simply speaking, the case of the tree we studied in this paper is an extrapolation of these two points depending on the conditions imposed at exits. When flows are imposed at outlets then each outlet is in Situation (i) and the optimal geometry is the one that tries to open the branches (in the limit of the constraint on the volume). When pressures are imposed, then all outlets are in Situation (ii) except one which is in Situation (i) since we imposed a non zero flow in the root of the tree and this flow, by conservation, has to get out of the tree by at least one branch. Moreover, Figure 12: Convergence curves of the Lagrange multiplier µ k , the Lagrangian L b (Ω k , µ k ), the viscous energy J(Ω k ) and the volume V (Ω k ). The x-coordinates represent the iterations number. The dashed line represents the step 4500 at which a remeshing was done.
in this last situation, our results show that it is "better" to close a branch and to use its volume to widen another branch. This implies that the other branch can distribute a larger amount of flow while dissipating less energy than two separate branches. Indeed, the viscous effects are larger near the walls and one wide branch has less walls than two smaller branches. When flow is imposed in a branch, this phenomena is compensated by the fact that reducing the branch radius increases the viscous effects by increasing the velocity gradients in the flow. Consequently, the optimal geometry is a compromise. These phenomena are probably true whenever the inertial effects are present or not, as our numerical simulations confirmed partially (Reynolds 100).
Finally, in term of modeling, our work shows that very different optimal structures can be obtained by boundary conditions adjustment (tree or pipe). For organs, we could imagine that depending on the organ function, optimization has been made through evolution by minimizing identical costs with different boundary conditions.

A Model proofs A.1 Proof of Proposition 1
The linearity of the relation between p and q comes from Poiseuille's law. Hence it is sufficient to compute the pressure vectors p related to the elements of canonical basis of R 2 N ×1 .
Let us begin with q = (1, 0, . . . , 0) ⊤ ∈ R 2 N ×1 . It corresponds to the case in which the fluid exits the tree only by the outlet denoted by (N, 1). In other terms, the fluid flows through the tree using the path Π 0→ (N,1) . By conservation, the volumetric flow rate which enters the tree through the root branch is exactly Φ = 1, so that the pressure p 1 at the outlet of the root branch is p 0 − r 0 . As there is no flow in the right-hand subtree stemming from the root branch, the pressure at its outlets (whose couples of indexes are between (N, 2 N −1 + 1) and (N, 2 N )) is the same as at the outlet of the root branch, namely p 0 − r 0 . Similarly, the pressure p 1,1 at the outlet of the branch denoted by (1, 1) is p 0 − (r 0 + r 1,1 ) = p 0 − r 0 1 + 1 ξ 1,1 . Following this approach recursively, one finds pressures at the outlet of the branches of the path Π 0→(N,1) denoted by (i, 1) with i ∈ 1, N to be p 0 − r 0 + (k,l)∈Π 0→(N,1) (i) r k,l , which writes again p 0 − r 0 1 + (k,l)∈Π 0→(N,1) (i) 1 ξ k,l . Therefore, the pressure p N,j (with j ∈ 1, 2 N ) at the outlet of the tree is Applying the same reasoning to any vector q = (0, . . . , 0, 1, 0, . . . , 0) ⊤ ∈ R 2 N ×1 (with 1 at the i-th position) yields The set of vectors (v i ) i∈ 1,2 N −1 forms a basis of {u N } ⊥ , the orthogonal complement of u N . Let i ∈ 1, 2 N − 1 . Multiplying (10) by v i in the sense of the inner product yields Furthermore, by conservation of the flow rate, one has q, u N = Φ. Thus, since A N (ξ) is real symmetric, the vector q verifies the following system q, A N (ξ)v i = − p, v i , ∀i ∈ 1, 2 N − 1 q, u N = Φ.
That explains the structure of the matrix M N (ξ). Let us now to prove that M N (ξ) is invertible. It amounts to prove that the set of vectors (A N (ξ)v 1 , . . . , A N (ξ)v 2 N −1 , u N ) is linearly independent. The fact that the set (A N (ξ)v 1 , . . . , A N (ξ)v 2 N −1 ) is linearly independent is almost trivial, since A N (ξ) is invertible by Proposition 2 and since the vectors (v i ) i∈ 1,2 N −1 are linearly independent. Now, assume that Multiplying the equation (53) by the vector A N (ξ) −1 u N in the sense of the inner product yields Since v i ∈ {u N } ⊥ for all i ∈ 1, 2 N − 1 , it follows that Moreover, u N , A N (ξ) −1 u N > 0. Indeed, we know that A N (ξ) is symmetric positive definite, and therefore A N (ξ) −1 inherits from these properties. Consequently, λ = 0. Now, (53) rewrites, By the linear independence of the set (A N (ξ)v i ) i∈ 1,2 N −1 , we conclude that ∀i ∈ 1, 2 N −1 , µ i = 0. It proves that the vectors (A N (ξ)v 1 , . . . , A N (ξ)v 2 N −1 , u N ) are linearly independent.

B An augmented Lagrangian algorithm
We present in this section the augmented Lagrangian algorithm used in this work to optimize the energy dissipated by a fluid in a dyadic tree with respect to the shape. One of the main interest of the augmented Lagrangian if compared to the classical Lagrangian is the regularization operation which generally improves the condition number of the dual function. Thus, a faster convergence of the maximizing sequence of the Lagrange multipliers is expected. However, the choice of a good parameter of augmentation can be very difficult. We will discuss this aspect concerning our simulations at the end of this section.
The descent direction in the main step of this algorithm will be computed thanks to a gradient method, which implies to calculate at each iteration the derivative of our criterion with respect to the domain.
The augmented Lagrangian associated to Problem (45) is where ℓ ∈ R is the Lagrange multiplier, associated with the volume constraint G(Ω), that is G(Ω) = meas(Ω) − V 0 .
We give now some precisions on the second step of the augmented Lagrangian algorithm, in particular on the choice of the descent method. Let Ω k be the domain obtained at iteration k−1, Γ k its lateral boundary and µ k the associated Lagrange multiplier. Ω k+1 is searched as a perturbation of the identity. That is why we write Ω k+1 = (I + ε k d k )(Ω k ), where d k is a vector field representing the perturbation of the mesh and ε k a variable step.
(d) Determination of the displacement of the mesh d k as the solution of the elliptic equation (e) Determination of an ε k that decreases the augmented Lagrangian.
3. Stopping criterion. The algorithm stops if |ℓ k+1 − ℓ k | ≤ ε stop and loops back to step 2 if the inequality is false.
As pointed out upwards, the choice of the parameter b in the augmented Lagrangian algorithm can be difficult. Practically speaking, b has to be chosen neither too big nor too small. Indeed, the biggest is the parameter b, the best is the conditioning of the dual functional giving the constraints and then the convergence of the sequence of Lagrange multipliers. Nevertheless, if b is chosen too big, the conditioning of the primal problem min{L b (Ω, ℓ k ), Ω ∈ E} deteriorates and it becomes more difficult to solve. Thus a compromise has to be done. Practically, a lot of preliminary tests have to be done to find b before the algorithm could be run properly.