PRESSURE BOUNDARY CONDITIONS FOR VISCOUS FLOWS IN THIN TUBE STRUCTURES: STOKES EQUATIONS WITH LOCALLY DISTRIBUTED BRINKMAN TERM

. The steady state Stokes-Brinkman equations in a thin tube structure is considered. The Brinkman term diﬀers from zero only in small balls near the ends of the tubes. The boundary conditions are: given pressure at the inﬂow and outﬂow of the tube structure and the no slip boundary condition on the lateral boundary. The complete asymptotic expansion of the problem is constructed. The error estimates are proved. The method of partial asymptotic dimension reduction is introduced for the Stokes-Brinkman equations and justiﬁed by an error estimate. This method approximates the main problem by a hybrid dimension problem for the Stokes-Brinkman equations in a reduced domain.

part of the thrombus behaves as a porous medium, but approaching the surface of the thrombus it corresponds better to the Newtonian fluid.
Also in the present paper we will consider the inflow-outflow boundary conditions involving pressure. These conditions are more convenient from the computational point of view. They also allow to compute the permeability of a piece of tissue containing a network of vessels. An extensive studying of the Stokes and Navier-Stokes equations with boundary conditions involving pressure started in the pioneer paper [8] and since then continues in a vast mathematical literature (see [35], [25], [17], [13], [16] and the bibliography there).

Definition of a thin tube structure
Let us remind the definitions of the tube structure and its graph given in [30]. Below in various situations, we choose one or another coordinates system denoting the local variable in both cases by x (e) and pointing out which end is taken as the origin of the coordinate system.
With every edge e j we associate a bounded domain σ j ⊂ R n−1 containing the origin and having C 2 -smooth boundary ∂σ j , j = 1, ..., M . For every edge e j = e and associated σ j = σ (e) we denote by Π (e) ε the cylinder Π (e) ε = x (e) ∈ R n : x (e) n ∈ (0, |e|), where x (e) = (x (e) 1 , ..., x (e) n−1 ), |e| is the length of the edge e and ε > 0 is a small parameter. Notice that the edges e j and Cartesian coordinates of nodes and vertices O j , as well as the domains σ j , do not depend on ε.
We will define as well a semi-infinite dilated cylinder Π x − O j ε ∈ ω j }. Every vertex O j is the end of one and only one edge e k which will also be denoted as e Oj ; we will re-denote as well the domain σ k associated to this edge as σ Oj . Notice that the subscript k may be different from j.  Definition 1.2. By a tube structure, we call the following domain Suppose that it is a connected set and that the boundary ∂B ε of B ε is C 2 -smooth except for the parts of the boundary which are the bases γ j ε = {x (e) ∈ σ Oj , x

Formulation of the problem. Existence and uniqueness of a solution
Let Γ = ∂B ε \ N j=N1+1 γ j ε be the lateral surface of the domain B ε . In the tube structure B ε we define the spaces W 1,2 (B ε ) = η ∈ W 1,2 (B ε ) : η| Γ = 0, η τ | γ j ε = 0, j = N 1 + 1, . . . , N , Let us consider the following boundary value problem for the steady-state Stokes equations in a tube structure divu = 0 in B ε , u = 0 on ∂B ε \ ∪ N j=N1+1 γ j ε , u τ = 0 on γ j ε , −ν∂ n u · n + p = c j /ε 2 on γ j ε , j = N 1 + 1, ..., N, where ν ε is the function (effective dynamical viscosity related to the porosity) equal to the positive constant ν (0) everywhere except for the balls B(O l , rε) and equal to the given functions ν ε (x) = ν (l) x−O l ε in the balls B(O l , rε), l = 1, ..., N 1 , K ε is the n × n symmetric matrix-valued function (inverse to effective permeability of porous medium) equal to zero everywhere except for the balls B(O l , rε) and equal to the given functions , ν ε is greater than some positive constant independent of ε, K ε ≥ 0 (non-negative matrix) ; ν (l) , K (l) are independent of ε; n is the unit normal vector to γ j ε , u τ = u − (u · n)n is the tangential component of the vector u, ∂ n g = ∇g · n is the normal derivative of g, c j are some given constants. This model was rigorously derived from the Navier-Stokes equation in a porous medium in [2] and was extensively studied in fluid mechanics [14]. This model describes the Newtonian flow in the tubes combined with the fluid filtration process through the zones B(O l , rε), simulating the eventual clots or thrombi. In these zones ν ε stands for the effective dynamical velocity taking into account the porosity of the clot, while K ε stands for the inverse to the effective permeability of the clot. The coefficients K ε , ν ε is supposed to be C 1 -smooth function, while ν ε is continuous in the closure of B ε .
In this section we prove the existence and uniqueness of the solution to problem (2.2) with the right-hand side f ∈ L 2 (B ε ). From the boundary condition u τ | γ j ε = 0 and the divergence equation divu = 0, it follows that −ν∂ n u · n| γ j ε = 0.Thus we can rewrite (2.2) with the right-hand side in the following form where p j stand for the constants c j /ε 2 , D = ∇ + ∇ T . Let us define a weak solution of problem (2.3) as a vector field u ∈ K 1,2 (B ε ) satisfying the integral identity for every η ∈ K 1,2 (B ε ). Here and below for any two n × n matrices A and B having entries a ij and b ij denote Consider an equivalent weak formulation: find a vector field u ∈ K 1,2 (B ε ) satisfying the integral identity for every η ∈ K 1,2 (B ε ). The equivalence of these formulations follows from the identity which is a corollary of the relation N j=N1+1 γ j ε η n dx = 0 for the solenoidal vector-valued function η. Let us explain this weak formulation heuristically; the rigorous analysis of the equivalence of the weak formulation and the classical one needs to study the regularity of the weak solution, see [9], [1] for the methods.
Identity (2.4) follows from equations (2.3) after multiplying them by η ∈ K 1,2 (B ε ) and integrating by parts in B ε . On the other hand, for a sufficiently regular solution u satisfying (2.4) there exists a pressure field p such that the pair (u, p) satisfies equations (2.3) 1,2 a.e. in B ε . Boundary conditions (2.3) 3,4,5 are satisfied in the sense of traces (see the definition of the space K 1,2 (B ε )). More exactly, function p is defined up to an additive constant but this constant can be chosen so that p satisfies (2.3) 5 . Indeed, take in (2.4) a smooth solenoidal function η satisfying the boundary conditions η| Γ = 0, η τ | γ j ε = 0, j = N 1 + 1, ..., N. Integrating by parts in Hence, there exists a pressure function p such that Let us fix arbitrary j ∈ {N 1 + 1, . . . , N }. Taking η ∈ K 1,2 0 (B ε ) such that η| γ k ε = 0 for k = j, we get The normal traces of functions η from K 1,2 0 (B ε ) satisfying the conditions η| γ k ε = 0 for k = j compile the whole space where c j is a constant (similarly to [16], [17]). Using these relations and taking now in (2.8) a test function Thus, Since the pressure p in the weak formulation is defined up to an additive constant, we may set c j = c N = 0, j = N 1 + 1, ..., N . Then from (2.9) we have These considerations justify the definition of the weak solution.
Proof. Define in K 1,2 (B ε ) the inner product [u, η] = Bε 1 2 ν ε (x)Du · Dη + K ε (x)u · η dx. Due to the Poincaré-Friedrichs inequality and the Korn inequality the corresponding norm is equivalent to the Dirichlet norm with a constant independent of ε. Using the Cauchy-Schwarz, the trace theorem and the Poincaré-Friedrichs inequality we derive the estimates (2.11) Finally, (2.12) Consider the linear functional Φ : Due to estimates (2.11) and (2.12) it is bounded in the Hilbert space K 1,2 (B ε ) and so, by the Riesz theorem there exists a unique function u ∈ K 1,2 (B ε ) such that for all η ∈ K 1,2 (B ε ), So, this function u is a unique solution to the weak formulation (2.3). Due to estimates (2.11) and (2.12) it satisfies estimate (2.10).
Remark 2.2. Notice also that the weak solution u of problem (2.3) belongs to the space W 2,2 (B ε ) whenever f ∈ L 2 (B ε ). The corresponding pressure belongs to W 1,2 (B ε ). This can be proved extending the solutions and the data by reflection over the sections γ j ε to a larger domain (see [9], [1]). Remark 2.3. We will also consider problem

13)
f i (x) = 0 in the neighborhood of bases γ j ε . A weak solution is defined as u ∈ K 1,2 (B ε ) satisfying for every η ∈ K 1,2 (B ε ). In this case Theorem 2.1 can be easily generalized with the estimate with the constant c independent of ε.
Let us now recover and estimate the pressure p corresponding to the weak solution u ∈ W 1,2 (B ε ) of problem (2.3). To do this, we need to prove that a linear bounded functional L defined on the space W 1,2 (B ε ) and vanishing on the subspace K 1,2 (B ε ) can be represented in the form . First of all, let us study the divergence equation in the space W 1,2 (B ε ) with the right-hand side from L 2 (B ε ). First of all recall the well-known result on the divergence equation.
Let Ω be a bounded domain in R n , n = 2, 3. Denote Consider the following problem (divergence equation): for given h ∈ L 2 (Ω), find a vector field w ∈W 1,2 (Ω) satisfying the equation Its formulation follows [31] (Lem. 3.7). However, in [31] there is a misprint in the constant which we fixed here. Its proof is in Appendix B.
Let σ ⊂ R n−1 be a bounded domain with Lipschitz boundary and Π = x : x ∈ σ, 0 < x n < 1 , be a cylinder in R n . First, we consider the divergence equation in the cylinder Π.
Lemma 2.7. Let h ∈ L 2 (Π). Then the divergence equation (2.20) admits at least one solution w ∈ W 1,2 (Π) vanishing on the part of the boundary ∂Π \ {x : x n = 0} and w τ | xn=0 = 0. The solution satisfies the estimate Proof. Let us extend h as an odd function to the larger cylinder Π = x : x ∈ σ, −1 < x n < 1 (with respect to x n ). Consider in Π the following divergence equation:
Then, according to Lemma 2.4, there exists the solution W ∈W 1,2 (B ε ) satisfying the estimate Without loss of generality we can assume that W has the odd component W and even component W n (with respect to x n ). Indeed, if we have some solution W to this divergence equation (2.22), we can consider the function W defined by relations Then W is still a solution to the divergence equation (2.22) and it satisfies the above parity conditions with respect to x n . This implies that W = 0 on the section x n = 0, that is W τ = 0 for x n = 0.
The restriction w of W on Π is a solution to problem (2.20).
Lemma 2.8. Let h be a function from L 2 (B ε ). Then the divergence equation admits at least one solution w ∈ W 1,2 (B ε ), satisfying the estimate Proof. Let us fix a cylinder Π (e) ε having the base γ N ε . Let us cut this cylinder at the distance ε from the base According to Lemma 2.6, it satisfies the estimate Now, in the cylinder C ε we construct a solution w 0 ∈ W 1,2 (C ε ) of the divergence equation satisfying w τ | γ N ε = 0 and the estimate The first inequality in (2.29) follows from Lemma 2.7 contracting 1/ε-times the cylinder Π in (2.21), the second inequality follows from (2.27). Since, by the construction of Lemma 2.7, the vector-function w 0 vanishes on the second base of the cylinder C ε , it can be extended by zero into B ε \ C ε . So, w 0 ∈ W 1,2 (B ε ). Finally, taking w = w 1 + w 0 , we finalize the proof.
Proof. In the proof of Lemma 2.6 for any h ∈ L 2 (B ε ) there was constructed a function w, belonging the space and estimate (2.19). One can easily check that this construction defines a bounded linear operator M −1 from Let us define the equivalence classes in W 1,2 (B ε ). We say that two functions w and v belong to the same class if they have the same divergence: divw = divv. We will call functions w belonging to the class W as representatives of this class. Define the sum W + V of two classes W and V as the equivalence class containing the function w + v, where w ∈ W and v ∈ V are the representatives of W and V respectively. Also, define the product αW of the equivalence class W by a real number α as the equivalence class containing the function αw, where w ∈ W is a representative of the class W . So, we can consider the vector space of the equivalence classes (known in the literature as the quotient space). Introduce the inner product in this space: if W and V are two classes and w ∈ W , v ∈ V , then the inner product in this quotient space is defined as (2.31) One can easily check that this definition satisfies the axioms of the inner product (bilinearity, symmetry and positivity of the associated norm) and is stable with respect to the choice of the representatives w and v of the classes W and V respectively.
is uniquely defined for all functions of the class W , and so, one can consider Φ as a linear functional on the vector space of equivalence classes. As the functional Φ is bounded with respect to the norm of W 1,2 (B ε ) and the operator M −1 is also bounded, where c ε is a constant depending on ε. So, Φ is bounded on the space of equivalence classes with inner product (2.31) and so, according to the Riesz theorem, there exists a unique equivalence class U such that for a representative u of U and for any w ∈ W 1,2 (B ε ) it can be represented in a form of an inner product (2.32) Taking now p = divu we complete the proof of the existence of p. Its uniqueness follows from the uniqueness of the equivalence class U , so that for all u from U divu is the same. Now we can introduce another weak formulation for problem (2.3), namely, formulation "with pressure": find u ∈ W 1,2 (B ε ) and p ∈ L 2 (B ε ) such that holds for every η ∈ W 1,2 (B ε ).
Proof. Applying Theorem 2.9 to the functional , we get the existence and uniqueness of the pressure p ∈ L 2 (B ε ) such that So, u and p satisfy integral identity (2.33). Evidently, applying estimates (2.11), (2.12), we get: Using Lemma 2.8, we can take in (2.33) the test function η such that divη = p(x) and Then (2.36) yields: and from (2.10) we get (2.34).
Note that u is the same in both weak formulations.
In this case a weak solution is defined as u ∈ K 1,2 (B ε ) satisfying for every η ∈ K 1,2 (B ε ). Now, Theorem 2.1 can be easily generalized with the estimate with the constant c independent of ε.

Asymptotic expansion of the solution
In this section we describe the construction of the asymptotic expansion. Let ζ ∈ C 2 (R) be even function independent of ε such that, ζ(t) = 0 if |t| ≤ 1/3, and ζ(t) = 1 if |t| ≥ 2/3. Denote e = e Oj (the edge with the end O j ) and x (e) the Cartesian coordinates corresponding to the origin O j and the edge e, i.e., is the orthogonal matrix relating the global coordinates x with the local ones x (e) .
The asymptotic expansion of the velocity field is sought in the form: where the first sum is taken over all edges having a vertex as an end point (and with the origin of the local coordinate system at the vertex), the second sum is taken over all remaining edges, and all terms in these sums are extended by zero out of cylinders Π  (they will be defined below), they are expanded in powers of ε; functions U [BLO l ,J] are the boundary layer type functions exponentially decaying at the infinity, they as well are expanded in powers of ε: The asymptotic expansion of the pressure for every half-cylinder Π (e) ε , x n < |e|/2, corresponding to the edge e = O l O i l , l = N 1 + 1, ..., N , (O l is the origin of the local coordinate system) is sought in the form: and on every half-bundle HB O l , l = 1, ..., N 1 , (O l is the origin of the local coordinate system) we define: n + a (e) − a (es) + a (es) where the terms of the sum are extended by zero out of cylinders Π and Here e s is the selected edge of the bundle (arbitrary chosen among edges of the bundle) and the local coordinates 0 (y (e) )dy (e) . Solve the conductivity problem on the graph for the function p 0 : (0) ) t (y (e) ). (2.48) Introduce the notation: where functions ν (l) and K (l) are extended out of the ball B(0, R) by ν (0) and zero respectively. Now to compensate the residual due to the multiplication of the Poiseuille flow by a cut-off function, we consider the boundary layer correctors and problems for these correctors set in a dilated bundle of semi-infinite cylinders Ω l : For l = 1, ..., N 1 the boundary layer problem for (U (y)) is: .
Consider now the conductivity problem of rank 1 on the graph for the function p 1 : 1 between the edges e and e s of the bundle. This problem as well has a unique solution p 1 . Note that problem (2.54) can be reduced to a problem on the graph with a right-hand side for a continuous unknown function as in [32]. Now, constants s Suppose that all terms of expansion (2.41)-(2.46) corresponding to the rank less or equal to j − 1 are known, and that the macroscopic pressure p j on the graph is known as well. Let us describe the passage from the rank j − 1 to the rank j.
Step 1. As the macroscopic pressure on the graph p j is known, define for every edge e constants s (2.55) Step 2. The boundary layer solution is a pair U solving the following Stokes system in Ω l , l = 1, ..., N 1 : j (y (e) , t) . in the sense of (2.53); these constants may be different for different outlets. Since the pressure function is defined up to an additive constant, we can fix the limit constant equal to zero for the outlet corresponding to the selected edge e s .
Step 3. Solve the conductivity problem on the graph for the function p Step 4. Finally, we find the pressure P [BLO l ] j (y) in boundary layer problem (2.56), (2.57) for l = 1, ..., N 1 : This step finalizes the passage from j to j + 1.

Residual
Consider the asymptotic expansion u (J) , p (J) of order J (see (2.41), (2.43)). By construction Here the first line of the right-hand side is the residual of the term K (l) u (J) , the second line of the righthand side comes from the pressure gradient term; this term is the only one that was not compensated by the boundary layer-in-space problems. The last line is the residual generated by the multiplication of the boundary layer correctors by the cut-off function ζ |x−O l | |e|min . Notice that terms appearing in this last line exponentially vanish because in the set supp 1 − ζ |x−O l | |e|min (where χ = 0) the order of this term in L 2 -norm is O(e −c1/ε ) with some positive constant c 1 (see the Appendix in [33]). From the obtained formulas it follows that hold.
holds for every η ∈ W 1,2 (B ε ). Applying Theorem 2.10 with estimate (2.34) we get: Evaluating now the norm of the difference u (J) and u (J+2) we obtain: Replacing J by J + 2 in (2.66) yields: So, the triangle inequality gives the first estimate (2.65). Applying the same argument to p (J) and p (J+4) , we get the second estimate (2.65).
Remark 2.13. The asymptotic expansion (2.41)-(2.46) can be slightly modified without loss of the accuracy.
Namely, the argument |x − O l | |e| min of the cut-off function ζ may be replaced by where the constant C J is chosen in such a way that Here Ω l,R = Ω l ∪ {|y| > R}. Indeed, by Theorem A1 [33], U [BLO l ,J] (y) and P [BLO l ,J] (y) exponentially vanish as |y| → ∞. Thus, So, the estimate is true for δ such that i.e. for δ = C J ε| ln ε|.
Notice that the constant C J =c 1 J +c 2 , wherec 1 ,c 2 are constants independent of J. Then the W 1,2 (B ε )-norm of the difference of the constructed asymptotic expansion u (J) and the modified one is of order O(ε J+(n−1)/2 ).

Method of asymptotic partial decomposition of the domain for the inflow-outflow boundary condition involving pressure
Using the obtained results we introduce and justify the method of asymptotic partial decomposition of the domain (MAPDD) for problem (2.2). This method first published in [27] reduces the dimension on the main part of B ε replacing the solution by the Poiseuille type flow and keeps the original full dimension in small neighbourhoods of the nodes and vertices. By this, it reduces the computational costs and accelerates the traditional numerical strategies.
Let us describe the algorithm of the MAPDD for the Stokes problem set in a tube structure B ε . Let δ be a small positive number much greater than ε (it will be chosen of the order O(ε| ln ε|)). For any edge e = O i O j of the graph (O i , O j are two nodes) introduce two hyperplanes orthogonal to this edge and crossing it at the distance δ from its ends. If in the edge e = O i O j only one end O i is a node and the second end is a vertex, then we introduce only one hyperplane at the distance δ from the node O i . Denote the cross-sections of the cylinder Π (e) ε by these two hyperplanes respectively, by S i,j (the cross-section at the distance δ from O i ) and S j,i (the cross-section at the distance δ from O j ) and denote the part of the cylinder Π Define the subspace W dec (B ε , δ) of W 1,2 (B ε ) such that on every truncated cylinder B dec,ε ij its elements (vector-valued functions) have on B dec,ε ij vanishing tangential entries of the vector-function and independent of the normal variable normal component of the vector-function. Let H dec (B ε , δ) be the subspace of W dec (B ε , δ) consisting of the divergence free functions. We will consider as well the subspace L 2 dec (B ε , δ) of the space L 2 (B ε ) such that its elements are affine functions of x (e) 1 on every truncated cylinder B dec,ε ij . The MAPDD approximation to problem (2.2) is formulated as a projection of the weak formulation (2.4) with f = 0 on the subspace H dec (B ε , δ), namely find v d ∈ H dec (B ε , δ), such that ∀η ∈ H dec (B ε , δ), the following integral identity holds.
Applying the standard Lax-Milgram lemma arguments one can prove the existence and uniqueness of the solution v d to this problem.
Note that the corrected, according to the above Remark 1.1, asymptotic solution u (J) belongs to the space H dec (B ε , δ) and still satisfies the weak formulation (2.4) with the same residual as before, i.e., of order O(ε J−2 ) in L 2 (B ε )-norm: Of course, this identity is still true for η ∈ H dec (B ε , δ), because H dec (B ε , δ) is a subspace of K 1,2 (B ε ). Evidently, the difference v d − u (J+2) belongs to H dec (B ε , δ) and for every η ∈ H dec (B ε , δ) satisfies the integral identity This estimate justifies the MAPDD.
(2.74) Also, introduce the space H dec ((U B) ε,δ ) as the subspace of divergence free functions of the space W dec ((U B) ε,δ ). Then problem (2.69) can be stated as follows: find v d ∈ H dec ((U B) ε,δ ), such that ∀η ∈ H dec ((U B) ε,δ ), the following integral identity holds. Here d l = |e| − 2δ, the distance between the cross-sections S ij and S ji cutting the segment e (for segments with the both end points which are nodes, e connects two nodes), d l = |e| − δ if one of the end points is a vertex, i.e. e connects a node and a vertex. In the last case, in (2.74) with j = N 1 + 1, ..., N or l = N 1 + 1, ..., N , by convention S lj = γ j ε (respectively, S lj = γ l ε ). The advantage of this form of the integral identity is that it uses only functions defined in small truncated domains.
It is possible to "recover" the MAPDD pressure and get an equivalent formulation "with pressure": find the vector-field v d and the "MAPDD pressure" p d for all i = 1, ..., N 1 belonging to the spaces v d ∈ W dec ((U B) ε,δ ) and p d ∈ L 2 (B ε,δ i ), and satisfying for all η ∈ W dec ((U B) ε,δ ) the following integral identity (2.76) Theorem 2.16. There exists a unique solution to problem (2.76) such that the vector field v d is the same as in Theorem 2.14, and the uniquely defined in all domains B ε,δ i pressure p d can be extended on all cylinders B dec,ε ij as an affine function of the local variable x (e) n so that the following error estimate holds. Here δ = C J+2n+1 ε| ln ε| with constant C J+2n+1 satisfying (2.68) with J replaced by J + 2n + 1.
The details of the proof of this theorem are presented in the Appendix A.

Some comments on the mixed boundary conditions
Consider now the case when on some part of the surfaces γ j ε , j = N 1 + 1, ..., N 2 , one sets the pressure conditions (2.2) 5 while on the other part (j = N 2 + 1, ..., N ) the inflow/outflow velocity is given. This case is treated in the same way as above. The boundary layers are constructed as above. There is no need of the compatibility condition on the integrals of the boundary value function g · n over γ j ε , j = N 2 + 1, ..., N (if N 2 < N 1 ). The equations on the graph have then boundary conditions for fluxes at the vertices O j , j = N 2 + 1, ..., N , and the boundary conditions for the unknown macroscopic pressure p at the vertices O j , j = N 1 + 1, ..., N 2 . Also note, that the results can be generalized to the case when ν ε and K ε are less regular: boundad measurable functions.

On the Darcy law for a tissue with network of vessels
The constructed above asymptotic expansion of the solution to problem (2.2) can be applied to the determining the permeability of a piece of tissue containing a network of vessels. The derivation of the Darcy law for flows in porous media with periodic structure from the Stokes and Navier-Stokes equations was introduced in [22], [10], [39]. Its justification was based on the method developed in [38], [40], [41] and the Appendix by L.Tartar in the book [39] (see also [22] for non-stationary setting and [43] and [4] for asymptotic expansions of the solution).
Consider a domain G containing a tube structure B ε . For simplicity, assume that G is a cube (0, 1) n and that all vertices and the surfaces γ j ε , j = N 1 + 1, ..., N belong to the faces of the cube S 0 = {x 1 = 0, (x 2 , x 3 ) ∈ (0, 1) n−1 } and S 1 = {x 1 = 1, (x 2 , x 3 ) ∈ (0, 1) n−1 }. Let all constants c j corresponding to γ j ε ⊂ S 1 are equal to ε 2 , so that p j = 1, while all remaining constants c j = 0. The cube G can be considered as a porous medium and according to the Darcy's law (confirmed by the above asymptotic expansion) the pressure gap (here equal to one) is proportional to the average velocity in the direction x 1 , equal toū 1 = Bε u 1 (x)dx. According to the above asymptotic analysis, this integral is of order ε n+1 because the velocity magnitude is O(ε 2 ) and the measure of the domain B ε is of order ε n−1 . Using the leading term of the asymptotic expansion ε 2Ũ (e) 0 (x (e) /ε) in each cylinder Π (e) ε and replacing the integration over B ε by the integration over M j=1 Π (ej ) ε we modify the average velocity with an error of order ε n+2 because the measure of all smoothing domains ω ε j is of order ε n , and the integral of the exponentially decaying boundary layer correctors is of order ε n+2 . Therefore, the algorithm to compute the permeability in x 1 direction is as follows: solve the problem on the graph (2.47), compute the approximate average velocityū where cos(n, x 1 ) is the cosine between the edge e j and the axis x 1 , and define the permeability in x 1 direction as the average velocity divided by the pressure gap (equal to one), so, the permeability is equal toū a 1 . Note that forū a 1 we get the following expression: Here s

Conclusion
The stationary Stokes equation with Brinkman term in small parts of the domain is studied in thin tube structures with the pressure condition at the inflows and outflows and no-slip boundary condition on the lateral boundary. This boundary value problem models the blood flow in a network of blood vessels, where the Brinkman flow zones simulate the filtration of the blood through thrombi. Also, it can be used to model the flow through a roll of thin capillaries, part of the network. The leading term of the asymptotic expansion can be used to determine the permeability of the tissue. The obtained MAPDD model of the flow can be coupled with the diffusion-convection equations modeling the transport of cells or substances by the blood in the same way as it was done in [7]. Such coupling can be done as well for the non-Newtonian flows in a network of vessels.
Appendix A. Proof of the recovery of the pressure in the MAPDD problem with constant C * independent of ε and δ.  u 0 (y (e) )dy (e) and u 0 is a unique solution of the problem −ν (0) ∆ y (e) u 0 (y (e) ) = 1, y (e) ∈ σ (e) , u 0 (y (e) )| ∂σ (e) = 0, y (e) ∈ ∂σ (e) . where Ψ l are given real numbers. The existence and uniqueness of the solution to this problem is proved as in [32].
Relation between κ i (P (e) ) * 0, ..., 0, u 0,ε (x (e) ) . Note that the order of the measure of cross section of the cylinder is ε n−1 and the magnitude of u 0,ε (x (e) ) in the Poiseuille velocity is of order ε 2 (and respectively, ε 4 is the square of the magnitude of u 0,ε ) while its derivatives are of order ε. So, we have: Then we extend U i inside the domains B ε,δ k as an arbitrary function from W 1,2 (B ε,δ k ) with the given boundary values (we keep the same notation for the extended function). In particular, we can do it just by multiplication of the Poiseuille velocities U i (x (e) ) by cut-off functions η ε = 1 − ζ depending only on the longitudinal variable x 1 − δ ε |) and we obtain the same estimates in B ε,δ k , namely In fact, the last estimates will contain an extra factor ε but it doesn't improve the overall result: let us calculate, for example, the norm Here we used the estimates η ε = O(1), η ε = O(ε −1 ), so that η ε 2 L 2 ((0,ε)) = O(ε −1 ) and η ε 2 L 2 ((0,ε)) = O(ε). Summing up all these estimates for U i we get estimate (A.1). The proof of the Lemma is completed.
i , belonging to L 2 (B ε,δ i ) for all i = 1, ..., N 1 . Then there exists a vector-valued function U ∈ W dec (B ε , δ) such that There holds the estimate with some positive constant C independent of ε and δ. For δ = O(ε| ln ε|) we get Proof. Consider the sum where U i are functions constructed in the previous lemma. By Stokes formula, and so, Thus, Now we need to solve the divergence equation in B ε,δ k , i.e. to construct a function Θ ∈W 1,2 (B ε,δ k ) such that The existence of a solution of (A.9) follows from Lemma 2.7. However, in order to obtain an appropriate estimate (with a constant independent of ε and δ), we need to reduce this problem to the same problem in the δ −1 -times dilated domain This domain is a thin tube structure with ε/δ as a small parameter replacing former small parameter ε (recall that δ = cε| ln ε|). So, we can apply Lemma 2.7 and construct the solution satisfying the estimate in original variables: (A. 10) Here and below c is a generic constant independent of small parameters. Let us evaluate the norm Ψ W 1,2 (B ε,δ k ) . It is majorated by the sum (see the previous lemma). So, finally, (A.11) Let us take U = Ψ + Θ, where Θ is extended by zero to the cylinders B dec,ε ij . Then U = Ψ for the remaining part of the tube structure. Recall that divΨ = 0 for x ∈ B ε \ N1 i=1 B ε,δ i . So, using estimate (A.11) we finalize the proof. Now we are in position to prove the existence of the solution (v d , p d ) to problem (2.76) with the test functions from the space W dec (B ε , δ).
Introduce the space L 2 dec (B ε ) as the space of scalar functions q ∈ L 2 dec (B ε ), equal to zero on all cylinders B dec,ε ij . Theorem A.3. There exists a unique function p d ∈ L 2 dec (B ε ) satisfying integral identity (2.76).
Proof. The proof of this theorem repeats the proof of Theorem 2.9 where the spaces L 2 (B ε ), K 1,2 (B ε ), and W 1,2 (B ε ) are replaced by L 2 dec (B ε ), H dec (B ε , δ), and W dec (B ε , δ) respectively, while the bounded linear operator holds, and there exists an extension p d ∈ L 2 dec (B ε ) of p d such that the following estimate holds.
Consider the difference q = p (J ) − p d . Applying Lemma A.2, we can construct a function U ∈ W dec (B ε , δ) such that and with some constant C independent of ε, and δ. Taking U as a test function in (A.14) and in (2.76) with J instead of J, consider the difference between these integral identities. Denoting u = u (J ) − v d , we get Applying estimate (2.67) for p − p (J ) , and then the triangle inequality we prove estimate (A.12). Let us extend now p d to the cylinders B dec,ε ij as an affine function by the following formula: where for any function q ∈ L 2 (B ε ) we denote by q Sij the mean value of q in the cylinder C ij = {x ε , i = 1, . . . , M, γ M +k = ω k ε , k = 1, . . . , N 1 , and γ il = γ i ∩ γ l . Notice that mes γ il is of order ε n , n = 2, 3. For each pair of intersecting domains γ i and γ l we introduce a C 2 -regular function η il (x) defined on B ε such that supp η il ⊂ γ il , |η il (x)| ≤ c ε n , |∇η il (x)| ≤ c ε n+1 and of the tube structure B ε , i.e., c il = 1, if i = l and γ il is nonempty, and c il = 0 in remaining cases. The existence of the matrix {K il } i,l=1,...,M +N1 is proved by induction on the number M + N 1 . If M + N 1 = 2, then c 12 = c 21 = 1 and by virtue of (B.2), one can set K 12 = −K 21 = µ 1 , K 11 = K 22 = 0. Assume that the assertion is valid for M + N 1 = m − 1. We prove it for M + N 1 = m. Let us set K il = 0 for those i and l for which c il = 0. As it is well known (see [15]), from a connected graph with the number of vertices greater than two, one can always discard one vertex so that the remaining graph is connected. From here there follows the existence of an index λ, 1 ≤ λ ≤ m, such that the matrix {c il }, i = λ, l = λ, is the adjacency matrix of some connected graph with (m − 1) vertices. Without loss of generality one can assume that λ = m. We consider the m-th equation from (B.4): Since the initial graph is connected, there exists an index l 1 , such that c ml1 = 1. We set K ml1 = −K l1m = µ m , while K ml = −K lm = 0 for l = l 1 . The remaining equations from (B.4) will be written in the form: Indeed, In addition, as it is mentioned above, the matrix {c il } i,l=1,...,m−1 satisfies the induction hypothesis for M + N 1 = m − 1. Consequently, system (B.6) has a solution with the required properties. Note that the matrix {c il } i.l=1,...,M +N1 is independent of ε, and it can be proved by induction (as before) that |K il | ≤ c max |µ 1 (t)|, . . . , |µ M +N1 | with the constant c independent of ε. Therefore, |K il | ≤ cε (n−1)/2 h L 2 (Bε) .
Define on B ε the function Using the skew-symmetry of {K il } i,l=1,...,M +N1 , it is easy to see that supp η i ⊂ γ i and The functions h (i) have the following properties: (i) h (i) have the same regularity as h,  (2.16) in ω k ε with the right-sides h (k) (x). For w (ei) and w (k) hold estimates of Lemmas 2.5 and 2.4 (for the contracted 1/ε times domain Ω) , respectively. Extend the functions w (ei) and w (k) by zero into B ε and put According to (B.7), So, Therefore, using Lemmas 2.4 and 2.5 we obtain required estimate for w.