Time periodic Navier–Stokes equations in a thin tube structure

The time periodic Navier–Stokes equations are considered in the three-dimensional and two-dimensional settings with Dirichlet boundary conditions in thin tube structures. These structures are finite union of thin cylinders (thin rectangles in the case of dimension two), where the small parameter ε is the ratio of the hight and the diameter of the cylinders. We consider the case of finite or big coefficient before the time derivative. This setting is motivated by hemodynamical applications. Theorems of existence and uniqueness of a solution are proved. Complete asymptotic expansion of a solution is constructed and justified by estimates of the difference of the exact solution and truncated series of the expansion in norms taking into account the first and second derivatives with respect to the space variables and the first derivative in time. The method of asymptotic partial decomposition of the domain is justified for the time periodic problem.


Introduction
The present paper is motivated by the problem of computer modeling of the blood vessel network. Such modeling is important for optimization of the choice of serging strategy in the case of cardio-vascular or cerebrovascular problems. The existing one-dimensional models and codes cannot give the required accuracy in the neighborhoods of clot formation zones, stents, and bifurcation of vessels. On the other hand, the completely threedimensional computations are currently very time consuming and can be applied only for small parts of the blood circulation system. That is why we suggest the hybrid dimension models, combining the one-dimensional reduction in the regular zones with threedimensional zooms in small zones of singular behavior. This approach was addressed to the non-stationary initial boundary value problems for the Navier-Stokes equations in thin tube structures [17]. However, for the hemodynamical modeling, more natural are time-periodic settings. That is why we consider here the periodic in time problem for the Navier-Stokes equation, prove the existence and uniqueness theorems for this setting, and construct the asymptotic expansion of a solution with respect to the small parameter ε equal to the ratio of the diameter of vessels to their length. We consider the Navier-Stokes equations in a class of special domains called tube structures. These domains are connected finite unions of thin finite cylinders (thin rectangles in the 2D case respectively). Each such tube structure may be schematically represented by its graph: letting the thickness of tubes tend to zero, we find out that tubes degenerate to segments. We consider the problem in two different scalings: one of them is the same as in [17], while the other generates a big coefficient of order ε -2 of the time derivative of the velocity. The constructed asymptotic expansion is then used for the construction of a special numerical method combining one-dimensional description with three-dimensional zooms, the method of asymptotic partial decomposition of the domain (MAPDD). This method reduces the full geometry setting to the computations in neighborhoods of bifurcation zone of diameter of order ε| ln ε| as in [4] but in the time periodic setting. An alternative approach was developed by A. Quarteroni's team [5], but this other method of junction of one-dimensional and three-dimensional zones is different because it is based on multi-physics modeling: the one-dimensional hyperbolic equations and three-dimensional models are derived independently of conservation laws, and then coupling is based on the ideas of consistency of numerical schemes implementing these models. On the contrary, the MAPDD starts from the Navier-Stokes equation written everywhere in the blood flow area, it rigorously derives the one-dimensional Poiseuille type equations in the main part of the domain with three-dimensional zoom in small parts near the bifurcations of vessels and clot formation zones. It prescribes mathematically justified size of the zoomed areas and asymptotically exact junction conditions. Numerous computational tests show that the multi-physics approach with hyperbolic one-dimensional models is more convenient for the description of thick vessels (arteries), while the MAPDD model works better for small vessels such as arterioles.
Let us describe now two different scalings of the Navier-Stokes equations in small and very small arterioles or capillaries. The experimental data depend on numerous factors: whether human or animal blood system is considered, if it is healthy or ill. So we take some averaged data from [10,11]. The characteristic time (period) is 1 second, while the characteristic velocity is about 0.5 × 10 -3 m/sec. Consider two scalings for the characteristic length and diameter of vessels: (1) the length is 10 -3 m and the characteristic diameter 10 -4 m, (2) the length is 10 -2 m and the characteristic diameter 10 -3 m (in both cases ε = 0.1). Let us make the change of the space variable X = 10 -3 x in case (1), and X = 10 -2 x in case (2). Consider case (1). Making the change of velocity v = 10 -3 V and the change of pressure p = 10 3 P and taking into account that the dynamic viscosity of the blood is about 4 × 10 -3 Pa sec and its density is 10 3 kg/m 3 , we obtain in new variables the Navier-Stokes equation with all coefficients of order one:
Consider now case (2). Making the change of velocity v = 10 -4 V and the change of pressure p = 10 1 P, we obtain in new variables the Navier-Stokes equation with all coefficients except for the time derivative term of order one, while the coefficient of the time derivative is 10 2 , i.e., ε -2 : 10 2 ∂V ∂t -4 X V + 5(V, ∇ X )V + ∇ X P = 0, ∇ · V = 0. That is why we consider below two different settings: with the factor ε -2 and without it.
In the first part of the paper we study the existence and uniqueness of the time periodic Navier-Stokes equations for both settings and derive the a priori estimates. The proof of the existence, uniqueness, and a priori estimates of a solution in the time periodic setting differs significantly from the proof of this trinity in the case of initial conditions given in [17]. Moreover, in the present paper we obtain more precise estimates by means of a technique using the Stokes operator extension (see [21]) and the base of its eigenfunctions. However, this technique needs the C 2 -smoothness of the boundary. That is why we modify slightly the definition of the tube structure adding smoothing domains near the vertices of the graph of the tube structure. This modification allows as well to improve the estimates for the Jth partial sums of an asymptotic expansion of solution constructed in Sect. 6 and to obtain better estimates than in [17]. The construction of the asymptotic expansion of the solution and the MAPDD implementation for setting (1) is similar to [17] but with the time periodicity condition instead of the initial condition.
On the contrary, case (2) corresponds to a different expansion of the solution, and the algorithm of its construction is given in Sect. 6. Section 8 is devoted to the justification of the asymptotic expansion.

Tube structures. Graphs
For the reader's convenience, let us remind the definitions of the tube structure and its graph given in [17], see also [12]. 3, and e 1 , e 2 , . . . , e M be M closed segments each connecting two of these points (i.e., each e j = O i j O k j , where i j , k j ∈ {1, . . . , N}, i j = k j ). All points O i are supposed to be the ends of some segments e j . The segments e j are called edges of the graph. A point O i is called a node if it is the common end of at least two edges and O i is called a vertex if it is the end of the only one edge. Any two edges e j and e i can intersect only at the common node. The set of vertices is supposed to be non-empty.
Denote by B = M j=1 e j the union of edges and assume that B is a connected set. The graph G is defined as the collection of nodes, vertices, and edges (see Fig. 1 In what follows, in the various situations we choose one or another coordinate 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 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 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 ε.
Let O 1 , . . . , O N 1 be nodes and O N 1 +1 , . . . , O N be vertices. Let ω 1 , . . . , ω N be bounded independent of ε domains in R n ; introduce the nodal domains ω Every vertex O j is the end of one and only one edge e k , which is re-denoted as e O j ; we re-denote as well the domain σ k associated with this edge as σ O j . Notice that the subscript k may be different from j. Definition 1. 2 We call the following domain a tube structure: Suppose that it is a connected set and that the boundary ∂B ε of B ε is C 2 -smooth.
In what follows we use the following notation. Let V be a Banach space. The norm of the element u in the function space V is denoted by u V . Vector-valued functions are denoted by bold letters, and the spaces of scalar and vector-valued functions are not distinguished in notation. We use the standard notations for Sobolev and Hölder spaces. Let T be a positive number. The notation V per means that elements of the space V are T-periodic functions, i.e., u(·, t) = u(·, t + T). Without loss of generality we may assume T = 2π . We will need two spaces of periodic functions: L 2 per (0, 2π) andW 1,2 per (0, 2π), which are supplied by the inner product of L 2 (0, 2π) andW 1,2 (0, 2π), respectively.

Formulation of the problem
Consider in the tube structure B ε the time periodic boundary value problem for the Navier-Stokes equations (1.1) Assume that the fluid velocity g at the boundary ∂B ε has the following structure: g = 0 everywhere on ∂B ε except for the set the corresponding dilated part of the boundary, g ∈ C [ J+1 2 ]+1 (0, 2π; W 3/2,2 (∂B ε )). Denote e = e O j (the edge with the end O j ) and x (e) is the Cartesian coordinates corresponding to the origin O j and the edge e, i.e., x (e) = P (e) (x -O j ), P (e) is the orthogonal matrix relating the global coordinates x with the local ones x (e) , σ where n is the unit outward (with respect to B ε ) normal vector to γ j ε , y (e) = x (e) ε ,ĝ j (y (e) , t) = g j ((P (e) ) * y (e) , t). Since g(x, t) is time periodic, F j (t) also must be time periodic. Moreover, since we will need the divergence-free extension of g, we assume the compatibility condition for the flow rates F j (t): Let g be the divergence-free time periodic extension of the boundary function g (which we denote by the same symbol g, g ∈ C [ J+1 2 ]+1 ([0, 2π]; W 2,2 (B ε ))) satisfying, for all t ∈ [0, 2π], the following asymptotic estimates: where the constant c is independent of ε. Below we construct the special extension g in the form of asymptotic representation of the solution g such that the discrepancy of this extension in equations (1.1) is small. But first we consider the following variations problem: find a vector-field v = u + g with div u = 0, u ∈ L ∞ per (0, 2π;W 1,2 (B ε )), u t ∈ L 2 per (0, 2π; L 2 (B ε )) satisfying the integral identity for every divergence-free vector-field η η η ∈W 1,2 (B ε ). Here g is an arbitrary extension satisfying (1.4) and f is an arbitrary function such that where c is independent of ε.

Auxiliary results
In this section we prove some multiplicative inequalities in the tube structure B ε . First we construct two coverings of the domain B ε . Take domains A n ∈ (0, 2ε)}, j = N 1 + 1, . . . , N (i.e., when O j are vertices), and A (j) x (e k j ) n ∈ (0, 2ε)}, j = 1, . . . , N 1 (i.e., when O j are nodes), where the union over k j is taken over all edges of the bundle B (j) associated with the node O j .
In parallel with the covering we take the covering

Lemma 2.1 (Poincaré inequality) The following inequality
holds with the constant c independent of ε.
with the constant c independent of ε.
Proof In any bounded Lipschitz domain Ω the following inequality holds (see [9]): and, by the Young inequality, By scaling, it is easy to see that in any domain A ε from the covering A ε we get the estimate Summing the above inequalities over all domains A ε from the covering A ε , we get Putting now δ 2 = ( u 2 L 2 (B ε ) + ε 2 ∇u 2 L 2 (B ε ) ) u -2 L 2 (B ε ) yields (2.7).
Proof In any bounded Lipschitz domain Ω the interpolation inequality holds (see [9]): By scaling in ε-domain A ε (i.e., A ε in any direction is of "size" ε), we derive Summing the above inequalities over all domains A ε from the covering A ε , we get with the constant c independent of ε.
Proof By the multiplicative inequality in a bounded Lipschitz domain Ω (see [9]) and, by the Young inequality, Then, by scaling, we get, for any domains A ε from the covering A ε , the estimates and summing them over all A ε from A ε , we derive Taking in the last inequality gives (2.10).
with the constant c independent of ε. In particular, Proof In any bounded Lipschitz domains Ω the following multiplicative inequality holds (see Lemma 13.2 in [1]): By scaling it is easy to see that in A ε we have Taking the supremum over all A ε ∈ A ε , we obtain (2.11). Inequality (2.12) follows from (2.11), (2.1), and (2.9).
with the constant c independent of ε. In particular, (2.14) Proof In any bounded Lipschitz two-dimensional domain Ω the following interpolation inequality holds (see Lemma 13.2 in [1]): By scaling, we have Taking the supremum over all A ε ∈ A ε , we get (2.13).

Stokes operator
Consider in B ε the Dirichlet problem for the Stokes system and hence the estimate with the constant c independent of ε.
Making the change of variables x = ε -1 y, we transform A ε and A ε into the fixed (independent of ε) domains A and A. The Stokes problem in coordinates y takes the form ADN local estimates for elliptic problems (see [2]) yield the inequality (see [8]). Multiplying (3.4) by w and integrating by parts yields Therefore, From (3.5), using (3.6) and the Poincaré inequality, we derive Returning to coordinates x, we obtain Estimating the last term in the right-hand side of (3.8), by (3.2), we derive (3.3).
Problem (3.1) can be rewritten in the operator form (without loss of generality we sup- , and P is the Leray projection onto divergence-free vector fields. By the same notation we denote the Friedrichs extension of this operator to the whole space H(B ε ). is called the Stokes operator. It is known that (see [7,21]): (i) The Stokes operator has a discrete spectrum: . Then, by density arguments, it follows that Thus, Here, V * means the dual space to V . From (3.3) we get the estimate which together with (3.10) gives and with the constant c independent of ε.
Proof We prove the solvability of problem (1.5) by Galerkin method (see [7,21]). The main purpose is to obtain suitable a priori estimates. The remaining part is standard. If u is a weak solution, then taking in (1.5) η η η = u, we obtain Using (1.4) and the Poincaré inequality (2.1), we derive the estimate and hence, multiplying this relation by 2ε β and using the Poincaré inequality, we get Multiplying this inequality with omitted second term in the left-hand side by e c * ε β-2 t and integrating over t, we obtain Consider Galerkin approximations of the solution to problem (1.5) defined by the following system of ordinary differential equations: . From estimate (4.4), which remains valid for Galerkin approximations, it follows that, for every N , the map M : The continuity of M follows from the general theory of systems of ordinary differential equations. It can be proved as well directly using the arguments of the proof of uniqueness Theorem 4.2. Since M is continuous, this ensures the existence of a 2π -periodic solution to the Galerkin approximations (4.5) for each fixed N . Now we derive a set of a priori estimates for Galerkin approximations u (N) . Integrate (4.3) with respect to t. Using the periodicity Because of the Poincaré inequality and the mean value theorem for Lebesgue integrals, there exists a point t * ∈ [0, 2π] such that Without loss of generality we may assume that t * = 0 (if not, we can consider problem (1.5) on the interval [t * , t * + 2π] and reduce it, by change of variable t → tt * , to the interval [0, 2π]). Integrating (4.3) from 0 to t and using (4.7), we get Estimates (4.7) and (4.8) are valid for Galerkin approximations constructed using an arbitrary basis and for arbitrary bounded Lipschitz domains. In order to estimate the higher derivatives of u, we have to assume that ∂B ε ∈ C 2 , and as a basis we shall use the eigenfunctions of the Stokes operator.

Solvability of problem (1.5); the case n = 3
In this section we prove the existence of the unique weak solution of problem (1.5).
holds. If the constant in (1.6) is sufficiently small (independently of ε) in the case β = 0 or if β = 2, then there holds also the estimate Proof As in Theorem 4.1, we use Galerkin approximations. First, applying inequality (2.3) instead of (2.2), we prove, exactly in the same way as before, the existence of Galerkin approximations and the following estimate for them: In order to estimate the higher derivatives of u, we use as a basis the eigenfunctions of the Stokes operator. Taking in (4.5) Galerkin approximations with the basis {w k } ∞ k=1 , multiplying it by λ k γ (N) k (t), summing up the obtained equalities from k = 1 to k = N , and using properties of the Stokes operator, we obtain (here we again omit the subscript N ) Let us estimate the right-hand side of (5.4). Using (2.1), (2.4)-(2.10), (3.12), and (1.4) for n = 3, we obtain Substituting (5.5)-(5.8) into (5.4) and taking δ = 1 12 implies If ε is sufficiently small (ε 1/4 ≤ 1 4C 1 ), then and (5.9) yields Denoting z(t) = B ε |∇u| 2 dx, we rewrite (5.10) as or, equivalently, Integrating (5.11) by t and using (5.3), we obtain Here, as in Sect. 4, we assume that t * = 0. For sufficiently small ε in the case β = 2 or for sufficiently small constant in inequality (1.6) in the case β = 0, the following inequalities hold. Then (5.12) gives Let us estimate the norm of u (N) t . Taking in integral identity (1.5) η η η = u (N) t (more precisely, multiplying (4.5) by d dt γ (N) k and summing by k from 1 to N ), we obtain (omitting the subscript N ) Integrating (5.15) over [0,2π], using the periodicity condition, (1.4), and inequalities (3.12), (5.14), for sufficiently small ε, we derive Estimates (5.3), (5.14), and (5.16) ensure the convergence of the Galerkin approximations and guarantee the existence of the solution (see [7,21]).

Asymptotic expansion
Let us describe the procedure of constructing an asymptotic expansion of the solution to problem (1.1) in the case β = 2. The case β = 0 is completely similar to the asymptotic expansion constructed in [17] with only one difference that all functions depending on time are 2π -periodic (instead of being equal to zero in some neighborhood of t = 0).
First we solve the time-periodic problem on the graph and find the macroscopic pressure as a periodic in time function linear on every edge with respect to the longitudinal variable x (e) n . At the nodes it satisfies the Kirchhoff type junction conditions. This problem on the graph is the time-periodic analogue of the problem considered in [15]. This problem defines in every cylinder Π (e) ε the Poiseuille type velocity depending only on the transversal space variable x (e) of the tube. We multiply the Poiseuille type velocity and pressure in every cylinder by cut-off functions ζ equal to one in the main middle part of the cylinder and vanishing in some O(ε)-neighborhood of the nodes. This multiplication generates an important residual in the right-hand side of the Navier-Stokes equations, having a finite support belonging to a O(ε)-neighborhood of the nodes. Then we construct the boundary layer correctors, compensating this residual. These correctors are solutions to the Stokes equations in the dilated bundles of cylinders extended by outlets to infinity. In this context the asymptotic expansion of the velocity is constructed in the form (y, t) exponentially decay as |y| tends to infinity. The asymptotic expansion of the pressure has a similar form: The asymptotic solution is constructed by induction with respect to j. At the base (initial) step j = 0, we consider the following problem on the graph: find a function p 0 ∈ L 2 per (0, 2π; W 1,2 (B)) such that equations hold. Here, Ψ l (t) = γ l g l · n dS. Operator L (e) relates the pressure slope S and the flux H in an infinite cylindrical pipe with section σ (e) . Namely, consider the following periodic in time boundary value problem for the heat equation: for given S ∈ L 2 per (0, 2π), find V ∈ L 2 per (0, 2π;W 1,2 (σ (e) )) with ∂V ∂t ∈ L 2 per (0, 2π; L 2 (σ (e) )) such that ∂V ∂t y (e) , tν y (e) V y (e) , t = S(t), y (e) ∈ σ (e) , t > 0, V y (e) , t ∂σ (e) = 0, V y (e) , t = V y (e) , t + 2π and denote V y (e) , t dy (e) = H(t), L (e) is a bounded linear operator acting from L 2 per (0, 2π) to W 1,2 per (0, 2π) (see [3,6]). Denote MS = V. The existence of a solution to problem (6.3) is proved in [18]. Let us represent p (e) 0 in the form For every edge e i , define the Poiseuille type velocity V (e i ) 0 (y (e i ) , t) as a vector such that in the local coordinates its last (i.e., normal) component is Ms (e) 0 , while the tangential components are equal to zero.
Next, we find the boundary layer correctors (V ) as solutions of the periodic in time Stokes equations in the dilated domain: union of semi-infinite cylinders having the common node O l , and the corresponding ω l . Namely, let O l be a node which is the common end of edges e i 1 , . . . , e i m of the bundle B (l) . Define the semi-infinite cylinders Π + l,j s = y ∈ R n : P (e is ) y ∈ σ i s × (0, +∞) and the domain Ω l with m outlets to infinity corresponding to the node O l : We introduce the boundary layer pressure of the rank -1 as V (e) 0,n y (e) , t , y ∈ Ω l , where the local coordinates have the origin at O l andâ (e) 1 (t) is an unknown function. This problem is decomposed to two independent ones: first we solve it without the term containingâ (e) 1 (t) in the right-hand side for (V and find a solution V which tends to zero as |y| → ∞, whileP [BLO l ] 0 at each outlet Π + l,j tends to a constantâ l,j (t), except for the outlet corresponding to a selected edge e s where it tends to zero (this is possible because the pressure is defined up to an additive constant). Then we solve the following problem on the graph: find a function p (e) 1 ∈ L 2 per (0, 2π; W 1,2 (e)) such that equations Analogously, if O l is a vertex, the end of the edge e i , then we define the domain Ω l , corresponding to this vertex, as follows: Ω l = y ∈ R n : P (e i ) y ∈ σ i × (0, +∞) ∪ ω l , and the boundary layer problem has the form Because of condition (6.3) 3 , we have ) which exponentially tends to zero at infinity (see [17,19,20]).
Suppose that all terms of asymptotic expansion corresponding to the rank less or equal to j -1 are known and the pressure on the graph p j is known as well. Describe the passage from rank j -1 to the rank j.
Step 1. As pressure on the graph p j is known, define for every edge e functions s (e) j (t) and a This problem is solved in two steps: first we find the couple (V ) which is the solution of the same problem without the last term in the definition of f It has a unique (up to an additive constant in the pressure) solution V It is satisfied because p (e) j is a solution to the problem on the graph with this junction condition (by the inductive hypothesis). When V [BLO l ] exponentially tends to zero as |y| → ∞, the corresponding pressure functionP [BLO l ] j stabilizes in outlets at infinity to some constantsâ (e) l,j (t); 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 the selected edge e s . Define ϕ (e) l,j+1 (t) =â (e) l,j (t).

is a parameter) if and only if
Similarly, in every vertex O l , l = N 1 +1, . . . , N , we get for (V ) the Stokes problem in Ω l which is the same as in the case of nodes O l with only one difference: there is no summing over e : O l ∈ e in the right-hand sides of the equations.
Step 3. Solve the problem on the graph for the function p (e) j+1 , (j < J): The local coordinates x (e) are defined so that all of them have the same origin O l .
Step 4. Finally we find the pressure P For j = J the last sum is absent. The last step finalizes the passage from j to j + 1.

Justification of the asymptotics
Consider Navier-Stokes problem (1.1). As an extension of the boundary value g, we take the asymptotic approximation u (J) constructed in the previous sections and let p (J) be the corresponding asymptotic approximation for the pressure p. By construction u (J) satisfies conditions (1.4). Represent v, p as the sums v = u + u (J) , p = q + p (J) . Then u, u (J) ∈ L 2 per (0, 2π; W 2,2 (B ε )), u t , u (J) t ∈ L 2 per (0, 2π; L 2 (B ε )). The difference u = v -u (J) is divergence free, satisfies the periodicity condition, the boundary condition u(x, t)| ∂B ε = 0, and the integral identity for every η η η ∈ H(B ε ). The existence of the unique solution u of (8.1) follows from Theorems 4.1-5.2.
Remark 8.4 The asymptotic expansion (6.1)-(6.2) can be slightly modified without loss of accuracy. Namely, the argument |x-O l | |e| min in the cut-off function ζ may be replaced by |x-O l | δ , where δ = C J ε| ln ε||e| min and the constant C J will be chosen below.
Denote J = J + 2. Consider the boundary layer functions V [BLO l ,J ] and P [BLO l ,J ] . It follows that these functions F [BLO l ,J ] (F stands for V or P) and their derivatives decay exponentially as the space variable tends to infinity in the outlets. Thus, there exist positive constants c 1 , c 2 such that, for all t ∈ [0, 2π] and for sufficiently large R, the following inequality holds: where Ω R l = Ω l ∪ {|y| > R}. Therefore, if B l ε = {x ∈ B ε : |x -O l | ≥ C J ε| ln ε||e| min /3}, then making change of the variable y = x-O l ε in the above inequality and taking R = C J | ln ε||e| min /3, we get ≤ c 1 exp -c 2 C J | ln ε||e| min /3 = c 1 ε c 2 C J |e| min /3 .
The proof of this theorem repeats the beginning of the proof of Theorem 4.1 (derivation of estimates (4.6) and (4.8)), where the Galerkin approximations are constructed in the space H 1 0,div=0 (B ε , δ) instead of H 1 div=0 (B ε ).