# Steady-state navier–stokes equations in thin tube structure with the bernoulli pressure inflow boundary conditions: Asymptotic analysis

Rita Juodagalvyte, Grigory Panasenko, Konstantinas Pileckas

► To cite this version:

Rita Juodagalvyte, Grigory Panasenko, Konstantinas Pileckas. Steady-state navier–stokes equations in thin tube structure with the bernoulli pressure inflow boundary conditions: Asymptotic analysis. Mathematics, 2021, 9 (19), pp.2433. 10.3390/math9192433. hal-04134243

HAL Id: hal-04134243

https://hal.science/hal-04134243v1

Submitted on 10 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Steady-State Navier-Stokes Equations in Thin Tube Structure with the Bernoulli Pressure Inflow Boundary Conditions: Asymptotic Analysis

Rita Juodagalvytė<fnref n="1" />,<fnref n="2" />, Grigory Panasenko<fnref n="1" />,<fnref n="2" />,* and Konstantinas Pileckas<fnref n="1" />

**Abstract:** Steady-state Navier-Stokes equations in a thin tube structure with the Bernoulli pressure inflow-outflow boundary conditions and no-slip boundary conditions at the lateral boundary are considered. Applying the Leray-Schauder fixed point theorem, we prove the existence and uniqueness of a weak solution. An asymptotic approximation of a weak solution is constructed and justified by an error estimate.

**Keywords:** Navier-Stokes equations; Bernoulli pressure boundary condition; asymptotic approximation; quasi-Poiseuille flows; boundary layers

**MSC:** 35Q35; 76D07

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Juodagalvytė, R.; Panasenko, G.; Pileckas, K. Steady-State Navier-Stokes Equations in Thin Tube Structure with the Bernoulli Pressure Inflow Boundary Conditions: Asymptotic Analysis. Mathematics **2021**, *9*, 2433. https://doi.org/10.3390/math9192433

Academic Editor: Christopher Goodrich

Received: 30 July 2021
Accepted: 22 September 2021
Published: 30 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors. License MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The asymptotic behavior of solutions of partial differential equations in thin domains is extensively studied in the vast mathematical literature. In particular, the thin tube structures, introduced in [1], are considered as a geometrical model of a blood vessel network (see other approaches to the modeling of blood vessel networks [2,3]). For the steady-state Navier-Stokes equations in a network of thin tubes, an asymptotic expansion of the solution was firstly constructed in [1]. The small parameter was introduced as the ratio of the thickness to the length of tubes in the network. This asymptotic expansion was used to justify the method of asymptotic partial decomposition of the domain firstly introduced for the stationary Navier-Stokes equations in thin tube structures in [4]. This method allowed reducing the computational costs that the Navier-Stokes equations posed in thin tube structures. In particular, the full-dimensional computations are only needed in small neighborhoods of the junction of tubes, while in the largest part of the domain, the computations are one-dimensional. The non-stationary Navier-Stokes equations in such a domain were studied in [5]. However, in these papers, the inflows and outflows were described by the given velocity at the corresponding parts of the boundary. For numerical implementation, the boundary conditions involving pressure for outflow are more natural. That is why such conditions were extensively studied in mathematical literature (see [6,7]). In particular, ref. [7] studies the stationary Navier-Stokes equations in a tube structure (a bundle of three tubes) with the given pressure at the "free" ends of the tubes. It is well known that this formulation of the problem has a solution for small data only. Therefore, ref. [7] proves such theorem of existence and uniqueness and constructs a first-order asymptotic approximation.

In the present paper, we consider the inflow-outflow boundary conditions for the Bernoulli pressure. These non-linear boundary conditions were studied first in [8] for the non-stationary Navier-Stokes equations for a compressible fluid and in [6] for the

<footnote n="1">¹ Institute of Applied Mathematics, Vilnius University, Naugarduko Str. 24, 03225 Vilnius, Lithuania; rita.juodagalvyte@mif.vu.lt (R.J.); konstantinas.pileckas@mif.vu.lt (K.P.)</footnote>

<footnote n="2">² Institute Camille Jordan UMR (CNRS 5208), University Jean-Monnet, 23 Rue P. Michelon, 42023 Saint-Etienne, France</footnote>

<footnote>* Correspondence: grigory.panasenko@univ-st-etienne.fr</footnote>
stationary Navier–Stokes equations for an incompressible fluid with small data where the existence of the solution was proved. Since then, these conditions were considered in different contexts: in [9], these conditions were studied for arbitrary data in a finite pipe; in [10], a special model of the decomposition of the boundary value problem for the non-Newtonian flow with the Bernoulli boundary conditions and some special newly introduced interface conditions inside the domain was studied, and the existence of a weak solution to this problem was proved.

In this publication, we will construct an asymptotic expansion of a weak solution of the stationary Navier–Stokes equations in the whole thin tube structure with the Bernoulli boundary conditions for the inflows and outflows. As an auxiliary result, we prove the existence and uniqueness of the solution. It is proved by the same method as in [9] but taking into account the dependence of the domain on the small parameter. The dependence of an a priori estimate on the small parameter is addressed. The main difficulty is related to the construction of boundary layers corresponding to the non-linear boundary conditions and the justification of the asymptotic expansion.

The paper has the following structure. Section 1 is the Introduction. Section 2 recalls the main definitions for thin structures and formulates the Poincaré–Friedrichs and embedding inequalities in thin structures. The main problem is formulated in Section 3. Also, Section 3 is devoted to the mathematical analysis of the stationary Navier–Stokes equations in thin tube structures with Bernoulli’s pressure condition on the inflow–outflow parts of the boundary, while on the lateral part of the boundary, the no-slip boundary condition is set. Theorem 1 states the existence of a weak solution for an arbitrary right-hand side from $L^2$ and constant given Bernoulli’s pressure on the inflows and outflows. An a priori estimate is proved with a constant depending on $\epsilon$, and this dependence is studied. Theorem 2 states that for the data bounded by some constant, the solution is unique, and the solutions corresponding to two different sets of the data satisfy the stability estimate: the norm of the difference of the solutions is bounded by some constant multiplied by an appropriate norm of the difference of the right-hand sides.

In Section 4, an asymptotic expansion of the solution of the stationary Navier–Stokes equations in the thin structure is constructed. This construction uses the stabilization theorem of the Stokes equations in a cylinder with the no-slip conditions on the lateral boundary, with the proof provided in Appendix A. Compared to the Navier–Stokes problem with a given velocity on the whole boundary, the problem with Bernoulli’s conditions leads to the Dirichlet problem on the graph for the macroscopic pressure, while in the case of a given velocity, the problem on the graph is of Neumann’s type.

In Section 5, the residual is calculated and evaluated. Finally, in Section 6, we prove the error estimate for the asymptotic approximations. Section 7 is the Conclusion.

## 2. Thin Tube Structure

Let us remind the reader of the definitions of the tube structure and its graph given in [11].

**Definition 1.** Let $O_1, O_2, \dots, O_N$ be $N$ different points in $\mathbb{R}^n$, $n = 2, 3$, and $e_1, e_2, \dots, e_M$ be $M$ closed segments, each connecting two of these points (i.e., each $e_j = \overline{O_{i_j}O_{k_j}}$, where $i_j, k_j \in \{1, \dots, N\}, i_j \neq 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 only intersect at the common node. The set of vertices is supposed to be non-empty.

Denote $\mathcal{B} = \bigcup_{j=1}^{M} e_j$ the union of edges, and assume that $\mathcal{B}$ is a connected set (see Figure 1). The union of all edges that have the same end point $O_l$ is called the bundle $\mathcal{B}_l$. Figure 1 presents the graph as a union of edges $e_1, \dots, e_7$, where points $O_1, O_2, O_3$ are the nodes, and points $O_4, O_5, O_6, O_7, O_8$ are the vertices. Each point $O_i$, a node or a vertex, with
all edges containing $O_i$ as an end point, form bundle $B_1$, for example, $O_1$ with edges $e_1, e_2$ and $e_5$ form bundle $B_1$.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Graph of the tube structure.</figcaption></figure>

Let $e$ be some edge, $e = \overline{O_iO_j}$. Consider two Cartesian coordinate systems in $\mathbb{R}^n$. The first one has the origin in $O_i$, and the axis $O_i x_n^{(e)}$ has the direction of the ray $[O_iO_j)$; the second one has the origin in $O_j$, and the opposite direction, i.e., $O_j \tilde{x}_n^{(e)}$, is directed over the ray $[O_jO_i)$.

Below, in various situations, we choose one or another coordinate's 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 $\sigma^j \subset \mathbb{R}^{n-1}$ containing the origin $O_i$ and having $\mathbb{C}^2$-smooth boundary $\partial\sigma^j$, $j = 1, \dots, M$. For every edge $e_j = e$ and associated $\sigma^j = \sigma^{(e)}$, we denote the cylinder by $\Pi_{\varepsilon}^{(e)}$

$$ \Pi_{\varepsilon}^{(e)} = \{ x^{(e)} \in \mathbb{R}^n : x_n^{(e)} \in (0, |e|), \frac{x^{(e)\prime}}{\varepsilon} \in \sigma^{(e)} \}, $$

where $x^{(e)\prime} = (x_1^{(e)}, \dots, x_{n-1}^{(e)})$, $|e|$ is the length of the edge $e$, and $\varepsilon > 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 $\sigma^j$, do not depend on $\varepsilon$. We will also define a semi-infinite dilated cylinder $\Pi_{\infty}^{(e)} = \{ x^{(e)} \in \mathbb{R}^n : x_n^{(e)} \in [0, \infty), x^{(e)\prime} \in \sigma^{(e)} \}$.

Let $O_1, \dots, O_{N_1}$ be nodes and $O_{N_1+1}, \dots, O_N$ be vertices. Let $\omega^1, \dots, \omega^{N_1}$ be the bounded domains independent of $\varepsilon$ in $\mathbb{R}^n$; introduce the nodal domains $\omega_{\varepsilon}^j = \{ x \in \mathbb{R}^n : \frac{x - O_j}{\varepsilon} \in \omega^j \}$.

Every vertex $O_j$ is the end of one and only one edge $e_k$, which will also be denoted as $e_{O_j}$; we will also re-denote the domain $\sigma^k$ associated to this edge as $\sigma^{O_j}$. Notice that the subscript $k$ may be different from $j$.

**Definition 2.** By a tube structure, we call the following domain

$$ B_{\varepsilon} = \left( \bigcup_{j=1}^{M} \Pi_{\varepsilon}^{(e_j)} \right) \cup \left( \bigcup_{j=1}^{N_1} \omega_{\varepsilon}^{j} \right). $$

Suppose that it is a connected set and that the boundary $\partial B_{\varepsilon}$ of $B_{\varepsilon}$ is $\mathbb{C}^2$-smooth except for the parts of the boundary that are the bases $\gamma_{\varepsilon}^j = \{ x^{(e)\prime} \in \sigma^{O_j}, x_n^{(e)} = 0 \}$ of cylinders $\Pi_{\varepsilon}^{(e)}$, $j = N_1 + 1, \dots, N$ (see Figure 2).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Tube structures.</figcaption></figure>

Let $r_1$ be the maximal diameter of domains $\omega^i$, $i = 1, ..., N$, denote $r = r_1 + 1$. Consider a node or a vertex $O_l$ and all edges $e_j$ having $O_l$ as one of their end points. We call the union of all these edges a bundle of edges and denote it $B_l$, i.e., $B_l = \bigcup_{j:O_l \in e_j} e_j$. By a bundle of cylinders $B_{O_l}$, we call the union $\omega^l_\epsilon \cup (\bigcup_{j:O_l \in e_j} \Pi^{(e_j)}_\epsilon)$. We will also consider the half-bundle $HB_{O_l} = \omega^l_\epsilon \cup (\bigcup_{j:O_l \in e_j} \{x \in \Pi^{(e_j)}_\epsilon, x_n^{(e_j)} \in [0, |e_j|/2]\})$. We will also use $\Omega_l = \omega^l \cup (\bigcup_{j:O_l \in e_j} \Pi^{(e_j)}_\infty)$, a bundle of dilated cylinders. Denote also $\Omega^\epsilon_l = \{x \in \mathbb{R}^n : x/\epsilon \in \Omega_l\}$.

Below, we will 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 $\Gamma = \partial B_\epsilon \setminus \bigcup_{j=N_1+1}^N \gamma_j^i$ be the lateral surface of the domain $B_\epsilon$. We introduce the function spaces $\hat{W}_\gamma^{1,2}(B_\epsilon)$ and $\hat{J}_\gamma^{1,2}(B_\epsilon)$:

$$ \begin{aligned} \hat{W}_{\gamma}^{1,2}(B_{\epsilon}) &= \{\eta \in W^{1,2}(B_{\epsilon}) : \eta|_{\Gamma} = 0, \eta|_{\gamma_j^i} = 0, j = N_1 + 1, \dots, N\}, \\ \hat{J}_{\gamma}^{1,2}(B_{\epsilon}) &= \{\eta \in \hat{W}_{\gamma}^{1,2}(B_{\epsilon}) : \mathrm{div}\eta = 0\}. \end{aligned} $$

We also introduce a subspace $J_\gamma^{1,2}(B_\epsilon)$ of $\hat{J}_\gamma^{1,2}(B_\epsilon)$ defined by

$$ J_{\gamma}^{1,2}(B_{\epsilon}) = \{\eta \in \hat{J}_{\gamma}^{1,2}(B_{\epsilon}) : \int_{\gamma_{\epsilon}^j} \eta_{\mathbf{n}} dS = 0, j = N_1 + 1, \dots, N\}. $$

**Lemma 1** (Poincaré–Friedrichs inequality). *The following inequality*

$$ \|u\|_{L^2(B_\epsilon)} \le C_{PF}\epsilon \|\nabla u\|_{L^2(B_\epsilon)}, \quad \forall u \in \hat{W}^{1,2}(B_\epsilon) \quad (1) $$

holds with the constant $C_{PF}$ independent of $\epsilon$.

The proof follows from the standard Poincaré–Friedrichs inequality in a bounded domain with the Dirichlet boundary condition on a part of the boundary, the partition of $B_\epsilon$ with parts of a diameter of the order of $\epsilon$ and the scaling argument.
**Lemma 2.** Let $B_{\varepsilon} \subset \mathbb{R}^n$, $n = 2,3$, $u \in \hat{W}^{1,2}(B_{\varepsilon})$. Then

$$ \begin{aligned} \|u\|_{L^4(B_\varepsilon)}^4 &\le c\varepsilon^{-2}\|u\|_{L^2(B_\varepsilon)}^2 (\|u\|_{L^2(B_\varepsilon)}^2 + \varepsilon^2\|\nabla u\|_{L^2(B_\varepsilon)}^2) \\ &\le c\varepsilon^2\|\nabla u\|_{L^2(B_\varepsilon)}^4, \quad n=2 \end{aligned} \quad (2) $$

and

$$ \begin{aligned} \|u\|_{L^4(B_\varepsilon)}^4 &\le c\varepsilon^{-3}\|u\|_{L^2(B_\varepsilon)} (\|u\|_{L^2(B_\varepsilon)}^2 + \varepsilon^2\|\nabla u\|_{L^2(B_\varepsilon)}^2)^{3/2} \\ &\le c\varepsilon\|\nabla u\|_{L^2(B_\varepsilon)}^4, \quad n=3 \end{aligned} \quad (3) $$

with the constant c independent of $\varepsilon$.

This lemma is proved in [12].

**Lemma 3.** Let $B_{\varepsilon} \subset \mathbb{R}^n$, $n = 2,3$, $u \in \hat{W}^{1,2}(B_{\varepsilon})$. Then

$$ \|u\|_{L^2(\gamma_{\varepsilon}^j)}^2 \le c\varepsilon^{-1} (\|u\|_{L^2(B_{\varepsilon})}^2 + \varepsilon^2 \|\nabla u\|_{L^2(B_{\varepsilon})}^2) \le c\varepsilon \|\nabla u\|_{L^2(B_{\varepsilon})}^2 \quad (4) $$

with the constant c independent of $\varepsilon$.

**Proof.** The inequality in Equation (4) follows immediately from well-known trace estimate

$$ \|v\|_{L^2(\partial\Omega)} \le c\|v\|_{W^{1,2}(\Omega)} $$

and scaling argument. $\square$

**3. Formulation of the Problem. Existence, Uniqueness and Stability of a Solution**

Let us consider the following boundary value problem for the steady-state Navier-Stokes equations in a tube structure $B_{\varepsilon}$

$$ \left\{ \begin{array}{l} -\nu\Delta\mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{u} + \nabla p = 0 \quad \text{in } B_{\varepsilon}, \\ \operatorname{div}\mathbf{u} = 0 \quad \text{in } B_{\varepsilon}, \\ \mathbf{u} = 0 \quad \text{on } \partial B_{\varepsilon} \setminus \bigcup_{j=N_1+1}^N \gamma_{\varepsilon}^j, \\ \mathbf{u}_{\tau} = 0 \quad \text{on } \gamma_{\varepsilon}^j, \\ -\nu\partial_n\mathbf{u} \cdot \mathbf{n} + \left(p + \frac{1}{2}|\mathbf{u}|^2\right) = c^j/\varepsilon^2 \quad \text{on } \gamma_{\varepsilon}^j, \quad j = N_1+1, \dots, N, \end{array} \right. \quad (5) $$

where $\nu$ is a positive constant, $\mathbf{n}$ is the unit normal vector to $\gamma_{\varepsilon}^j$, $\mathbf{u}_{\tau} = \mathbf{u} - (\mathbf{u} \cdot \mathbf{n})\mathbf{n}$ is the tangential component of the vector $\mathbf{u}$, $\partial_n g = \nabla g \cdot \mathbf{n}$ is the normal derivative of $g$, $c^j$ are some constants.

Note that in the fifth equation the first term is equal to zero.

In this section, we prove the existence and uniqueness of the solution to Equation (5), with the right-hand side $\mathbf{f} \in L^2(B_{\varepsilon})$. From the boundary condition $\mathbf{u}_{\tau}|_{\gamma_{\varepsilon}^j} = 0$ and the divergence equation $\operatorname{div}\mathbf{u} = 0$, it follows that $-\nu\partial_n\mathbf{u} \cdot \mathbf{n}|_{\gamma_{\varepsilon}^j} = 0$. Thus, using the identity $\frac{1}{2}(\nabla\mathbf{u}^2) = \mathbf{u} \cdot (\nabla\mathbf{u})^t = \sum_{i=k}^n u_k \nabla u_k$, we can rewrite Equation (5) with the right-hand side in the following form

$$ \left\{ \begin{array}{l} -\nu\Delta\mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{u} - \mathbf{u} \cdot (\nabla\mathbf{u})^t + \nabla\Phi = \mathbf{f} \quad \text{in } B_{\varepsilon}, \\ \operatorname{div}\mathbf{u} = 0 \quad \text{in } B_{\varepsilon}, \\ \mathbf{u} = 0 \quad \text{on } \partial B_{\varepsilon} \setminus \bigcup_{j=N_1+1}^N \gamma_{\varepsilon}^j, \\ \mathbf{u}_{\tau} = 0 \quad \text{on } \gamma_{\varepsilon}^j, \\ \Phi = p_j \quad \text{on } \gamma_{\varepsilon}^j, \quad j = N_1+1, \dots, N, \end{array} \right. \quad (6) $$

where $\Phi = (p + \frac{1}{2}|\mathbf{u}|^2)$ is the Bernoulli pressure, $p_j$ stand for the constants $c^j/\varepsilon^2$.
Let us define a weak solution of problem (6) as a vector field $\mathbf{u} \in \hat{J}_{\gamma}^{1,2}(B_{\varepsilon})$ satisfying the integral identity

$$ \begin{aligned} & \nu \int_{B_{\varepsilon}} \nabla \mathbf{u} : \nabla \eta \, dx + \int_{B_{\varepsilon}} (\mathbf{u} \cdot \nabla) \mathbf{u} \cdot \eta \, dx - \int_{B_{\varepsilon}} (\eta \cdot \nabla) \mathbf{u} \cdot \mathbf{u} \, dx \\ & = - \sum_{j=N_1+1}^{N} p_j \int_{\gamma_{\varepsilon}^j} \eta_{\mathbf{n}} \, dx' + \int_{B_{\varepsilon}} \mathbf{f} \cdot \eta \, dx \end{aligned} \quad (7) $$

for every $\eta \in \hat{J}_{\gamma}^{1,2}(B_{\varepsilon})$.

Introduce $p_j^* = p_j - p_N$, $j = N_1, ..., N$. Consider an equivalent weak formulation: find a vector field $\mathbf{u} \in \hat{J}_{\gamma}^{1,2}(B_{\varepsilon})$ satisfying the integral identity

$$ \begin{aligned} & \nu \int_{B_{\varepsilon}} \nabla \mathbf{u} : \nabla \eta \, dx + \int_{B_{\varepsilon}} (\mathbf{u} \cdot \nabla) \mathbf{u} \cdot \eta \, dx - \int_{B_{\varepsilon}} (\eta \cdot \nabla) \mathbf{u} \cdot \mathbf{u} \, dx \\ & = - \sum_{j=N_1+1}^{N-1} p_j^* \int_{\gamma_{\varepsilon}^j} \eta_{\mathbf{n}} \, dx' + \int_{B_{\varepsilon}} \mathbf{f} \cdot \eta \, dx \end{aligned} \quad (8) $$

for every $\eta \in \hat{J}_{\gamma}^{1,2}(B_{\varepsilon})$. The equivalence of these formulations follows from the identity

$$ \sum_{j=N_1+1}^{N-1} p_j^* \int_{\gamma_{\varepsilon}^j} \eta_{\mathbf{n}} \, dx' = \sum_{j=N_1+1}^{N} p_j \int_{\gamma_{\varepsilon}^j} \eta_{\mathbf{n}} \, dx', $$

which is a corollary of the relation

$$ \sum_{j=N_1+1}^{N} p_j^* \int_{\gamma_{\varepsilon}^j} \eta_{\mathbf{n}} \, dx' = 0 $$

for the solenoidal vector-valued function $\eta$.

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 [13] for the methods.

Identity (7) follows from Equation (6) after multiplying them by $\eta \in \hat{J}_{\gamma}^{1,2}(B_{\varepsilon})$ and integrating by parts in $B_{\varepsilon}$. On the other hand, for a sufficiently regular solution $\mathbf{u}$ satisfying Equation (7), there exists a pressure field $p$ such that the pair $(\mathbf{u}, p)$ satisfies Equation (6)$_{1,2}$ i.e., in $B_{\varepsilon}$. The boundary conditions (Equation (6))$_{3,4,5}$ are satisfied in the sense of traces (see the definition of the space $\hat{J}_{\gamma}^{1,2}(B_{\varepsilon})$). More exactly, function $\Phi$ is defined up to an additive constant, but this constant can be chosen so that $\Phi$ satisfies Equation (6)$_{5}$. Indeed, in Equation (7), take a smooth solenoidal function $\eta$ satisfying the boundary conditions $\eta|_{\Gamma} = 0$, $\eta|_{\gamma_{\varepsilon}^j} = 0$, $j = N_1 + 1, ..., N$. Integrating by parts in Equation (7) yields

$$ \begin{aligned} & \int_{B_{\varepsilon}} (-\nu \Delta \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} - \mathbf{u} \cdot (\nabla \mathbf{u})^t - \mathbf{f}) \cdot \eta \, dx \\ & = -\nu \sum_{j=N_1+1}^{N} \int_{\gamma_{\varepsilon}^j} \partial_n \mathbf{u} \cdot \eta \, dS - \sum_{j=N_1+1}^{N} p_j \int_{\gamma_{\varepsilon}^j} \eta_{\mathbf{n}} \, dx' \\ & = - \sum_{j=N_1+1}^{N} p_j \int_{\gamma_{\varepsilon}^j} \eta_{\mathbf{n}} \, dx'. \end{aligned} \quad (9) $$

If $\eta \in J_0^{\infty}(B_{\varepsilon}) = \{\eta \in C_0^{\infty}(B_{\varepsilon}) : \text{div}\eta = 0\}$, then it follows from Equation (9) that

$$ \int_{B_{\varepsilon}} (L(\mathbf{u}) - \mathbf{f}) \cdot \eta \, dx = 0 \quad \forall \eta \in J_0^{\infty}(B_{\varepsilon}), $$
where

$$L(\mathbf{u}) = -\nu \Delta \mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{u} - \mathbf{u} \cdot (\nabla \mathbf{u})^t.$$

Hence, there exists a pressure function $\Phi$, such that (e.g., [14])

$$L(\mathbf{u}) + \nabla\Phi = \mathbf{f} \quad \text{a.e. in } B_{\varepsilon}.$$

Then

$$\int_{B_{\varepsilon}} (L(\mathbf{u}) - \mathbf{f}) \cdot \boldsymbol{\eta} \, dx = - \int_{B_{\varepsilon}} \nabla\Phi \cdot \boldsymbol{\eta} \, dx = - \sum_{j=N_1+1}^{N} \int_{\gamma_{\varepsilon}^j} \Phi \cdot \boldsymbol{\eta}_n \, dx'$$

for every solenoidal $\boldsymbol{\eta} \in \hat{J}_{\gamma}^{1,2}(B_{\varepsilon})$. Therefore,

$$\sum_{j=N_1+1}^{N} \int_{\gamma_{\varepsilon}^j} \Phi \cdot \boldsymbol{\eta}_n \, dx' = \sum_{j=N_1+1}^{N} p_j \int_{\gamma_{\varepsilon}^j} \boldsymbol{\eta}_n \, dx'$$

Thus,

$$\sum_{j=N_1+1}^{N} \int_{\gamma_{\varepsilon}^j} (\Phi - p_j) \cdot \boldsymbol{\eta}_n \, dx' = 0 \quad \forall \boldsymbol{\eta} \in \hat{J}_{\gamma}^{1,2}(B_{\varepsilon}). \qquad (10)$$

Let us fix arbitrary $j \in \{N_1 + 1, \dots, N\}$. Taking $\boldsymbol{\eta} \in \hat{J}_{\gamma}^{1,2}(B_{\varepsilon})$ such that $\boldsymbol{\eta}|_{\gamma_{\varepsilon}^k} = 0$ for $k \neq j$, we obtain

$$(\Phi - p_j)|_{\gamma_{\varepsilon}^j} = c_j,$$

where $c_j$ is a constant (see [9,15]). Using these relations and now taking a test function $\boldsymbol{\eta} \in \hat{J}_{\gamma}^{1,2}(B_{\varepsilon})$ in Equation (10) such that $\int_{\gamma_{\varepsilon}^k} \boldsymbol{\eta}_n \, dx' = 0$ for $k \neq j$ and $k \neq N$, $\int_{\gamma_{\varepsilon}^j} \boldsymbol{\eta}_n \, dx' = 1$ and $\int_{\gamma_{\varepsilon}^N} \boldsymbol{\eta}_n \, dx' = -1$, we obtain

$$\sum_{j=N_1+1}^{N} c_j \int_{\gamma_{\varepsilon}^j} \boldsymbol{\eta}_n \, dx' = c_j - c_N \Rightarrow c_j = c_N.$$

Thus,

$$c_j = c_N \quad \forall j = N_1 + 1, \dots, N. \qquad (11)$$

Since the Bernoulli pressure $\Phi$ in the weak formulation is defined up to an additive constant, we may set $c_j = c_N = 0$, $j = N_1 + 1, \dots, N$. Then from Equation (11), we have

$$\Phi|_{\gamma_{\varepsilon}^j} = p_j, \quad j = N_1 + 1, \dots, N.$$

**Theorem 1.** For arbitrary $\mathbf{f} \in L^2(B_{\varepsilon})$ and $p_j^* \in \mathbb{R}$, $j = N_1 + 1, \dots, N - 1$ problem (6) admits at least one weak solution $\mathbf{u} \in \hat{J}_{\gamma}^{1,2}(B_{\varepsilon})$. There holds the estimate

$$\|\nabla \mathbf{u}\|_{L^2(B_{\varepsilon})} \le c \left( \varepsilon^{n/2} \sum_{j=N_1+1}^{N-1} |p_j^*| + \varepsilon \|\mathbf{f}\|_{L^2(B_{\varepsilon})} \right) \qquad (12)$$

with the constant $c$ independent of $\varepsilon$.
**Proof.** Define in $\hat{\mathcal{J}}_{\gamma}^{1,2}(B_{\varepsilon})$ the inner product $[\mathbf{u}, \eta] = \int_{B_{\varepsilon}} \nabla \mathbf{u} : \nabla \eta dx$ corresponding to the Dirichlet norm. Using the Hölder inequality and Lemmas 1 and 2, we derive the estimates

$$ \begin{aligned} & \left| \int_{B_{\varepsilon}} (\eta \cdot \nabla) \mathbf{u} \cdot \mathbf{u} dx \right| + \left| \int_{B_{\varepsilon}} (\mathbf{u} \cdot \nabla) \mathbf{u} \cdot \eta dx \right| \\ & \le \left( \int_{B_{\varepsilon}} |\mathbf{u}|^4 dx \right)^{1/4} \left( \int_{B_{\varepsilon}} |\nabla \mathbf{u}|^2 dx \right)^{1/2} \left( \int_{B_{\varepsilon}} |\eta|^4 dx \right)^{1/4} \\ & \le c \varepsilon^{\alpha} \|\nabla \mathbf{u}\|_{L^2(B_{\varepsilon})}^2 \|\nabla \eta\|_{L^2(B_{\varepsilon})}, \end{aligned} $$

where $\alpha = 1$ for $n = 2$ and $\alpha = 1/2$ for $n = 3$. From Lemma 3, it follows that

$$ \begin{aligned} & \left| \sum_{j=N_1+1}^{N-1} p_j^* \int_{\gamma_\varepsilon^j} \eta_n dx' \right| \le \sum_{j=N_1+1}^{N-1} |p_j^*| \left( \int_{\gamma_\varepsilon^j} |\eta|^2 dx \right)^{1/2} |\gamma_\varepsilon^j|^{1/2} \\ & \le c \varepsilon^{n/2} \sum_{j=N_1+1}^{N-1} |p_j^*| \|\nabla \eta\|_{L^2(B_\varepsilon)}. \end{aligned} \quad (13) $$

Finally,

$$ \begin{aligned} & \left| \int_{B_{\varepsilon}} \mathbf{f} \cdot \eta dx \right| \le \left( \int_{B_{\varepsilon}} |\mathbf{f}|^2 dx \right)^{1/2} \left( \int_{B_{\varepsilon}} |\eta|^2 dx \right)^{1/2} \\ & \le c \varepsilon \|\mathbf{f}\|_{L^2(B_{\varepsilon})} \|\nabla \eta\|_{L^2(B_{\varepsilon})}. \end{aligned} \quad (14) $$

From above estimates and the Riesz theorem, it follows that the integral identity (8) is equivalent to the operator equation in the space $\hat{\mathcal{J}}_{\gamma}^{1,2}(B_{\varepsilon})$:

$$ \mathbf{u} = \mathcal{A} \mathbf{u}, \quad (15) $$

where the operator $\mathcal{A}$ is defined by

$$ \begin{aligned} [\mathcal{A}\mathbf{u}, \eta] = & \int_{B_{\varepsilon}} v^{-1} \left[ -(\mathbf{u} \cdot \nabla)\mathbf{u} \cdot \eta + (\eta \cdot \nabla)\mathbf{u} \cdot \mathbf{u} + \mathbf{f} \cdot \eta \right] dx \\ & -v^{-1} \sum_{j=N_1+1}^{N-1} p_j^* \int_{\gamma_{\varepsilon}^j} \eta_n dx' \quad \forall \eta \in \hat{\mathcal{J}}_{\gamma}^{1,2}(B_{\varepsilon}). \end{aligned} $$

Using the compactness of the embedding $W^{1,2}(B_{\varepsilon}) \hookrightarrow L^4(B_{\varepsilon})$, it is standard to show that the operator $\mathcal{A}: \hat{\mathcal{J}}_{\gamma}^{1,2}(B_{\varepsilon}) \mapsto \hat{\mathcal{J}}_{\gamma}^{1,2}(B_{\varepsilon})$ is compact (see [14]). Thus, the existence of at least one solution to Equation (15) will follow from the Leray–Schauder fixed point theorem if we show that all possible solutions $\mathbf{u}^{(\lambda)}$ of the equation

$$ \mathbf{u}^{(\lambda)} = \lambda \mathcal{A} \mathbf{u}^{(\lambda)}, \quad \lambda \in [0, 1] \quad (16) $$

are uniformly (with respect to $\lambda$) bounded.

A solution $\mathbf{u}^{(\lambda)}$ to Equation (16) satisfies the integral identity

$$ \begin{aligned} & v \int_{B_{\varepsilon}} \nabla \mathbf{u}^{(\lambda)} : \nabla \eta dx + \lambda \int_{B_{\varepsilon}} (\mathbf{u}^{(\lambda)} \cdot \nabla) \mathbf{u}^{(\lambda)} \cdot \eta dx - \lambda \int_{B_{\varepsilon}} (\eta \cdot \nabla) \mathbf{u}^{(\lambda)} \cdot \mathbf{u}^{(\lambda)} dx \\ & = -\lambda \sum_{j=N_1+1}^{N-1} p_j^* \int_{\gamma_{\varepsilon}^j} \eta_n dx' + \lambda \int_{B_{\varepsilon}} \mathbf{f} \cdot \eta dx \quad \forall \eta \in \hat{\mathcal{J}}_{\gamma}^{1,2}(B_{\varepsilon}). \end{aligned} \quad (17) $$

Taking $\eta = \mathbf{u}^{(\lambda)}$ in Equation (17), we obtain

$$ v \int_{B_{\varepsilon}} |\nabla \mathbf{u}^{(\lambda)}|^2 dx = -\lambda \sum_{j=N_1+1}^{N-1} p_j^* \int_{\gamma_{\varepsilon}^j} \mathbf{u}_n^{(\lambda)} dx' + \lambda \int_{B_{\varepsilon}} \mathbf{f} \cdot \mathbf{u}^{(\lambda)} dx. $$
Using Equations (13) and (14), we obtain

$$ \|\nabla \mathbf{u}^{(\lambda)}\|_{L^2(B_\varepsilon)}^2 \le c \left( \varepsilon^{n/2} \sum_{j=N_1+1}^{N-1} |p_j^*| + \varepsilon \|\mathbf{f}\|_{L^2(B_\varepsilon)} \right) \|\nabla \mathbf{u}^{(\lambda)}\|_{L^2(B_\varepsilon)} $$

Hence

$$ \|\nabla \mathbf{u}^{(\lambda)}\|_{L^2(B_\varepsilon)} \le c \left( \varepsilon^{n/2} \sum_{j=N_1+1}^{N-1} |p_j^*| + \varepsilon \|\mathbf{f}\|_{L^2(B_\varepsilon)} \right). $$

The constant c in the last inequality is independent of λ and ε. This finishes the proof of the theorem. □

Define α = 1 for n = 2 and α = 1/2 for n = 3.

**Theorem 2.** 1. There exists a positive constant c₀ independent of ε, such that if

$$ c_0 \varepsilon^{\alpha} \left( \varepsilon^{n/2} \sum_{j=N_1+1}^{N-1} |p_j^*| + \varepsilon \|\mathbf{f}\|_{L^2(B_\varepsilon)} \right) < \nu, \quad (18) $$

then the weak solution **u** ∈ $\hat{J}_\gamma^{1,2}(B_\varepsilon)$ of Equation (6) is unique.

2. Let {$p_{1j}^*$} and {$p_{2j}^*$} $j = N_1 + 1, \dots, N$ be two sets of real constants and $\mathbf{f}_1, \mathbf{f}_2$ be functions, $\mathbf{f}_i \in L^2(B_\varepsilon), i = 1, 2$, satisfying Equation (18), and let $\mathbf{u}_i \in \hat{J}_\gamma^{1,2}(B_\varepsilon)$ be weak solutions of problem (6) corresponding to {$p_{ij}^*$} and $\mathbf{f}_i, i = 1, 2$. Then there exists a constant C independent of ε such that

$$ \|\nabla \mathbf{u}_1 - \nabla \mathbf{u}_2\|_{L^2(B_\varepsilon)} \le \varepsilon C \|\mathbf{f}_1 - \mathbf{f}_2\|_{L^2(B_\varepsilon)} + c \varepsilon^{n/2} \sum_{j=N_1+1}^{N-1} |p_{1j}^* - p_{2j}^*|. \quad (19) $$

**Proof.** 1. Suppose that there exist two weak solutions $\mathbf{u}_1$ and $\mathbf{u}_2$ satisfying Equation (8). Subtracting identity (8) for $\mathbf{u}_2$ from the one for $\mathbf{u}_1$, we obtain

$$ \nu \int_{B_\varepsilon} \nabla \mathbf{w} : \nabla \eta \, dx + \int_{B_\varepsilon} [(\mathbf{w} \cdot \nabla)\mathbf{u}_1 \cdot \eta + (\mathbf{u}_2 \cdot \nabla)\mathbf{w} \cdot \eta] \, dx \\ - \int_{B_\varepsilon} [(\eta \cdot \nabla)\mathbf{u}_1 \cdot \mathbf{w} + (\eta \cdot \nabla)\mathbf{w} \cdot \mathbf{u}_2] \, dx = 0, \quad (20) $$

where $\mathbf{w} = \mathbf{u}_1 - \mathbf{u}_2$. Taking $\eta = \mathbf{w}$ in Equation (20), in virtue of Lemma 2 and Equation (12) for $\mathbf{u}_2$, we obtain

$$ \begin{aligned} \nu \int_{B_\varepsilon} |\nabla \mathbf{w}|^2 dx &= \int_{B_\varepsilon} (\mathbf{w} \cdot \nabla) \mathbf{w} \cdot \mathbf{u}_2 dx - \int_{B_\varepsilon} (\mathbf{u}_2 \cdot \nabla) \cdot \mathbf{w} dx \\ &\le 2 \|\nabla \mathbf{u}_2\|_{L^4(B_\varepsilon)} \|\nabla \mathbf{w}\|_{L^2(B_\varepsilon)} \|\mathbf{w}\|_{L^4(B_\varepsilon)} \\ &\le c_0 \varepsilon^\alpha \left( \varepsilon^{n/2} \sum_{j=N_1+1}^{N-1} |p_j^*| + \varepsilon \|\mathbf{f}\|_{L^2(B_\varepsilon)} \right) \|\nabla \mathbf{w}\|_{L^2(B_\varepsilon)}^2, \end{aligned} $$

where the constant c₀ is independent of ε. If condition (18) is valid, the last inequality yields

$$ \int_{B_\varepsilon} |\nabla \mathbf{w}|^2 dx = 0, $$

and, thus, $\mathbf{w} = 0$.
2. Subtracting identity (8) for $\mathbf{u}_2$ from the one for $\mathbf{u}_1$, we obtain

$$ 
\begin{aligned}
& v \int_{B_{\varepsilon}} \nabla \mathbf{w} : \nabla \boldsymbol{\eta} \, dx + \int_{B_{\varepsilon}} [(\mathbf{w} \cdot \nabla)\mathbf{u}_1 \cdot \boldsymbol{\eta} + (\mathbf{u}_2 \cdot \nabla)\mathbf{w} \cdot \boldsymbol{\eta}] \, dx \\
& - \int_{B_{\varepsilon}} [(\boldsymbol{\eta} \cdot \nabla)\mathbf{u}_1 \cdot \mathbf{w} + (\boldsymbol{\eta} \cdot \nabla)\mathbf{w} \cdot \mathbf{u}_2] \, dx \\
& = \int_{B_{\varepsilon}} (\mathbf{f}_2 - \mathbf{f}_1) \cdot \boldsymbol{\eta} \, dx - \sum_{j=N_1+1}^{N-1} (p_{1j}^* - p_{2j}^*) \int_{\gamma_{\varepsilon}^j} \boldsymbol{\eta}_n \, dx',
\end{aligned} 
\quad (21) $$

where $\mathbf{w} = \mathbf{u}_1 - \mathbf{u}_2$. Taking $\boldsymbol{\eta} = \mathbf{w}$ in Equation (21), in virtue of Lemma 2 and Equation (12) for $\mathbf{u}_2$, we obtain

$$ 
\begin{aligned}
v \int_{B_{\varepsilon}} |\nabla \mathbf{w}|^2 dx &= \int_{B_{\varepsilon}} (\mathbf{w} \cdot \nabla) \mathbf{w} \cdot \mathbf{u}_2 dx - \int_{B_{\varepsilon}} (\mathbf{u}_2 \cdot \nabla) \mathbf{w} \cdot \mathbf{w} dx + \int_{B_{\varepsilon}} (\mathbf{f}_2 - \mathbf{f}_1) \cdot \mathbf{w} dx \\
&\quad - \sum_{j=N_1+1}^{N-1} (p_{1j}^* - p_{2j}^*) \int_{\gamma_{\varepsilon}^j} \mathbf{w}_n \, dx' \\
&\le c \|\nabla \mathbf{w}\|_{L^2(B_{\varepsilon})}^2 \|\mathbf{u}_2\|_{L^4(B_{\varepsilon})} + \|\mathbf{f}_2 - \mathbf{f}_1\|_{L^2(B_{\varepsilon})} \|\mathbf{w}\|_{L^2(B_{\varepsilon})} \\
&\quad + c \varepsilon^{n/2} \sum_{j=N_1+1}^{N-1} |p_{1j}^* - p_{2j}^*| \|\nabla \mathbf{w}\|_{L^2(B_{\varepsilon})}^2 \\
&\le c_0 \varepsilon^{\alpha} \left( \varepsilon^{n/2} \sum_{j=N_1+1}^{N-1} |p_j^*| + \varepsilon \|\mathbf{f}_2\|_{L^2(B_{\varepsilon})} \right) \|\nabla \mathbf{w}\|_{L^2(B_{\varepsilon})} \\
&\quad + \varepsilon C_{PF} \|\nabla \mathbf{w}\|_{L^2(B_{\varepsilon})} \|\mathbf{f}_2 - \mathbf{f}_1\|_{L^2(B_{\varepsilon})} + c \varepsilon^{n/2} \sum_{j=N_1+1}^{N-1} |p_{1j}^* - p_{2j}^*| \|\nabla \mathbf{w}\|_{L^2(B_{\varepsilon})}^2,
\end{aligned}
$$

where $\varepsilon C_{PF}$ is the Poincaré-Friedrichs constant for domain $B_{\varepsilon}$ (constant $C_{PF}$ is independent of $\varepsilon$). If condition (18) is valid, the last inequality yields Equation (19). □

**Remark 1.** Notice also that the weak solution $\mathbf{u}$ of problem (6) belongs to the space $W^{2,2}(B_{\varepsilon})$ whenever $\mathbf{f} \in L^2(B_{\varepsilon})$. The corresponding pressure belongs to $W^{1,2}(B_{\varepsilon})$. This can be proved extending the solutions and the data by reflection over the sections $\gamma_{\varepsilon}^j$ to a larger domain (see [13]).

### 4. Asymptotic Expansion of the Solution

In this section, we describe the construction of the asymptotic expansion. Let $\zeta \in C^2(\mathbb{R})$ be an even function independent of $\varepsilon$ such that $\zeta(t) = 0$ if $|t| \le 1/3$ and $\zeta(t) = 1$ if $|t| \ge 2/3$. Denote $e = e_{O_j}$ (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., $x^{(e)} = P^{(e)}(x - O_j)$, and $P^{(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:

$$ 
\begin{aligned}
\mathbf{u}^{(J)}(x) &= \sum_{O_l, l=N_1+1, \dots, N, \varepsilon=\overline{O_l O_{i_l}}} \zeta\left(\frac{|e| - x_n^{(e)}}{3r\varepsilon}\right) \mathbf{U}^{[e, J]}\left(\frac{x^{(e)\prime}}{\varepsilon}\right) \\
&\quad + \sum_{e=\overline{O_a O_b}; a, b \le N_1} \zeta\left(\frac{x_n^{(e)}}{3r\varepsilon}\right) \zeta\left(\frac{|e| - x_n^{(e)}}{3r\varepsilon}\right) \mathbf{U}^{[e, J]}\left(\frac{x^{(e)\prime}}{\varepsilon}\right) \\
&\quad + \sum_{l=1}^{N} \left(1 - \zeta\left(\frac{|x - O_l|}{|e|_{\min}}\right)\right) \mathbf{U}^{[BLO_l, J]}\left(\frac{x - O_l}{\varepsilon}\right),
\end{aligned} 
\quad (22) $$

where the first sum is taken over all edges that have a vertex as an end point (and with the origin of the local coordinate system at the vertex), and the second sum is taken over all
remaining edges. All the terms in these sums are extended by zero out of cylinders $\Pi_{\varepsilon}^{(e)}$; the terms of the third sum are extended by zero out of the corresponding bundles:

$$ 
\left\{ 
\begin{aligned}
\mathbf{U}^{[e,J]} &= (P^{(e)})^t (0, \dots, 0, \tilde{U}^{[e,J]})^t, \\
\mathbf{U}_{(j)}^{(e)} &= (P^{(e)})^t (0, \dots, 0, \tilde{U}_j^{(e)})^t, \quad j = 0, 1, \dots, J, \\
\tilde{U}^{[e,J]}(y^{(e)})' &= \sum_{j=0}^{J} \varepsilon^j \tilde{U}_j^{(e)}(y^{(e)})^t, \\
\mathbf{U}^{[BLO_l, J]}(y) &= \sum_{j=0}^{J} \varepsilon^j \mathbf{U}_j^{[BLO_l]}(y).
\end{aligned}
\right. 
\qquad (23) 
$$

The asymptotic expansion of the pressure for every half-cylinder $\Pi_{\varepsilon}^{(e)}$, $x_n < |e|/2$, corresponding to the edge $e = \overline{O_l O_{l'}}$, $l = N_1 + 1, \dots, N$, ($O_l$ is the origin of the local coordinate system) is sought in the form:

$$ p^{(J)}(x) = -s^{(e)} x_n^{(e)} + a^{(e)} + \frac{1}{\varepsilon} \left(1 - \zeta\left(\frac{|x - O_l|}{|e|_{min}}\right)\right) P^{[BLO_l, J]}\left(\frac{x - O_l}{\varepsilon}\right), \quad (24) $$

and on every half-bundle $HB_{O_l}$, $l = 1, \dots, N_1$, ($O_l$ is the origin of the local coordinate system) we define:

$$ p^{(J)}(x) = \sum_{e \subset B_l} \zeta\left(\frac{x_n^{(e)}}{3r\varepsilon}\right) \left(-s^{(e)}x_n^{(e)} + a^{(e)} - a^{(e_s)}\right) + a^{(e_s)} \\ + \frac{1}{\varepsilon}\left(1 - \zeta\left(\frac{|x - O_l|}{|e|_{min}}\right)\right) P^{[BLO_l, J]}\left(\frac{x - O_l}{\varepsilon}\right), \quad (25) $$

where the terms of the sum are extended by zero out of cylinders $\Pi_{\varepsilon}^{(e)}$,

$$ s^{(e)} = \frac{1}{\varepsilon^2} \sum_{j=0}^{J} \varepsilon^j s_j^{(e)}, \quad a^{(e)} = \frac{1}{\varepsilon^2} \sum_{j=0}^{J} \varepsilon^j a_j^{(e)} \qquad (26) $$

and

$$ P^{[BLO_l, J]}(y) = \sum_{j=0}^{J} \varepsilon^j P_j^{[BLO_l, J]}(y). \qquad (27) $$

Here $e_s$ is the selected edge of the bundle (arbitrary chosen among edges of the bundle), and the local coordinates $x^{(e)}$ are redefined so that all of them have the same origin $O_l$.

The algorithm of successive determination of the terms in asymptotic Equations (22) and (24) is as follows.

**The base case.** Solve the conductivity problem on the graph for the function $p_0$:

$$ 
\left\{ 
\begin{aligned}
& -\kappa_e \frac{\partial^2 p_0^{(e)}}{\partial x_n^{(e)2}} (x_n^{(e)}) = 0, \quad x_n^{(e)} \in (0, |e|), \\
& -\sum_{e:O_l \in e} \kappa_e \frac{\partial p_0^{(e)}}{\partial x_n^{(e)}}(0) = 0, \quad l = 1, \dots, N_1, \\
& p_0^{(e)}(0) = c^l, \quad l = N_1 + 1, \dots, N, \\
& p_0^{(e)}(0) = p_0^{(e_s)}(0), \quad \forall e \subset B_l.
\end{aligned}
\right. 
\qquad (28) 
$$

Here, the local coordinates $x^{(e)}$ are redefined so that all of them have the same origin $O_l$. Thus, $p_0$ is a continuous function on the graph. Indeed, the last condition of this problem means that the values of the function $p_0$ for the values of local variables $x_n^{(e)} = 0$ associated to all edges $e$ of the bundle $B_l$ are the same. Note that by applying the same
Lax-Milgram lemma arguments as in the first part of [16], one can prove the existence and uniqueness of the solution of this problem.

Solving the above conductivity problem, we define the constants $s_0^{(e)}$ and $a_0^{(e)}$ for every edge $e$ such that

$$p_0^{(e)}(x^{(e)}) = -s_0^{(e)}x_n^{(e)} + a_0^{(e)}$$

and the velocity

$$\tilde{U}_{(0)}^{(e)}(y^{(e)}) = s_{(0)}^{(e)} U_{0}^{(e)}(y^{(e)}), \quad \mathbf{U}_{0}^{(e)}(y^{(e)}) = (\mathcal{P}^{(e)})^{t}(0, \dots, 0, \tilde{U}_{(0)}^{(e)})^{t}(y^{(e)}), \quad (29)$$

where $U_0^{(e)}(y^{(e)})$ is the solution to the Dirichlet problem

$$-\nu \Delta_{(y^{(e)})} U_0^{(e)}(y^{(e)}) = 1, \quad y^{(e)} \in \sigma^{(e)}; \\ U_0^{(e)}|_{\partial\sigma^{(e)}} = 0.$$

For $l = 1, \dots, N_1$, the boundary layer problem for $(\mathbf{U}_0^{[BLO_l]}(y), P_0^{[BLO_l]}(y))$ is:

$$-\nu \Delta_y \mathbf{U}_0^{[BLO_l]} + \nabla_y P_0^{[BLO_l]} = \mathbf{f}_0^{[REGO_l]} + \mathbf{f}_0^{[BLO_l]}, \quad y \in \Omega_l, \\ \mathrm{div}_y \mathbf{U}_0^{[BLO_l]} = h_0^{[REGO_l]}, \quad y \in \Omega_l, \\ \mathbf{U}_0^{[BLO_l]}|_{\partial\Omega_l} = 0, \quad (30)$$

where

$$\mathbf{f}_0^{[REGO_l]}(y) = -\sum_{e:O_l \in e} \left\{ s_0^{(e)} \left( -\nu \Delta_y \left( \zeta \left( \frac{y_n^{(e)}}{3r} \right) (\mathcal{P}^{(e)})^t (0, \dots, 0, U_0^{[e]}(y^{(e)}))^* \right) \right. \right. \\ \left. \left. - \nabla_y \left( \zeta \left( \frac{y_n^{(e)}}{3r} \right) y_n^{(e)} \right) \right) + (a_1^{(e)} - a_1^{(e_s)}) \nabla_y \left( \zeta \left( \frac{y_n^{(e)}}{3r} \right) \right) \right\}, \\ \mathbf{f}_0^{[BLO_l]}(y) = 0, \quad (32)$$

$$h_0^{[REGO_l]}(y) = \mathrm{div}_y \sum_{e:O_l \in e} \left\{ s_0^{(e)} \zeta \left( \frac{y_n^{(e)}}{3r} \right) (\mathcal{P}^{(e)})^t (0, \dots, 0, U_0^{[e]}(y^{(e)}))^t \right\}. \quad (33)$$

Here the sum $\sum_{e:O_l \in e}$ is taken over all edges $e$ that have ends in the node $O_l$, and the terms are extended by zero out of each cylinder $\Pi_e^{(e_j)}$. Here we have an unknown quantity in the right-hand side, the constant $a_1^{(e)} - a_1^{(e_s)}$ is unknown. Let us denote by $(\mathbf{U}_0^{[BLO_l]}(y), \hat{P}_0^{[BLO_l]}(y))$ the solution of Equation (30) without the last term $(a_1^{(e)} - a_1^{(e_s)}) \nabla_y \left( \zeta \left( \frac{y_n^{(e)}}{3r} \right) \right)$ in $\mathbf{f}_0^{[REGO_l]}(y)$ (since this term is in the gradient form, the solutions only differ by the pressure components). The right-hand sides of Equation (30) have compact supports. Therefore, according to results of Propositions 7.1 and 7.2 in [5], the pressure $\hat{P}_0^{[BLO_l]}(y)$ exponentially stabilizes in each outlet (corresponding to the edge $e$) to a constant, say $\hat{a}_0^{[BLO_l,e]}$, in the sense of integral estimates

$$\lim_{k \to +\infty} \int_{\{y_n^{(e)} \in (k, k+1)\} \cap \Omega_l} (\hat{P}_0^{[BLO_l]}(y) - \hat{a}_0^{[BLO_l,e]})^2 dy = 0. \quad (34)$$

For $l = N_1 + 1, \dots, N$, $(\mathbf{U}_0^{[BLO_l]}(y), P_0^{[BLO_l]}(y)) = (0, 0)$.
Consider now the conductivity problem of rank 1 on the graph for the function $p_1$:

$$ 
\left\{ 
\begin{array}{l} 
-\kappa_e \frac{\partial^2 p_1^{(e)}}{\partial x_n^{(e)2}}(x_n^{(e)}) = 0, \quad x_n^{(e)} \in (0, |e|), \\ 
-\sum_{e:O_l \in e} \kappa_e \frac{\partial p_1^{(e)}}{\partial x_n^{(e)}}(0) = 0, \quad l = 1, \dots, N_1, \\ 
p_1^{(e)}(x_n^{(e)} = 0) = 0, \quad l = N_1 + 1, \dots, N, \\ 
p_1^{(e)}(0) - p_1^{(e_s)}(0) = \hat{a}_0^{[BLO_l, e]}, \quad \forall e \subset B_l, e \neq e_s, 
\end{array} 
\right. 
\quad (35) 
$$

where $e_s$ is the selected edge of the bundle. Therefore, in this problem on the graph, the solution may be discontinuous at the nodes. Namely, at each node $O_l$, there are prescribed jumps of $p_1^{(e)}$ between the edges $e$ and $e_s$ of the bundle. This problem also has a unique solution $p_1$.

Now, constants $s_1^{(e)}$ and $a_1^{(e)}$ are known: $p_1^{(e)}(x_n^{(e)}) = -s_1^{(e)}x_n^{(e)} + a_1^{(e)}$, and we can completely determine the pressure in the boundary layer problem (30):

$$ P_0^{[BLO_l]}(y) = \hat{P}_0^{[BLO_l]}(y) - \sum_{e:O_l \in e, e \neq e_s} \zeta \left( \frac{y_n^{(e)}}{3r} \right) \hat{a}_0^{[BLO_l, e]}. $$

Suppose that all terms of Equations (22)–(27) corresponding to the rank less or equal to $j-1$ are known, and the macroscopic pressure on the graph $p_j$ is known as well. Let us describe the passage from rank $j-1$ to rank $j$.

**Step 1.** As the macroscopic pressure on the graph $p_j$ is known, define constants $s_j^{(e)}$ and $a_j^{(e)}$ for every edge $e$ such that

$$ p_j^{(e)}(x^{(e)}) = -s_j^{(e)}x_n^{(e)} + a_j^{(e)} $$

and

$$ 
\begin{aligned}
\tilde{U}_{(j)}^{(e)}(y^{(e)}) &= s_j^{(e)} U_0^{(e)}(y^{(e)}), \\
\mathbf{U}_j^{(e)} &= (\mathcal{P}^{(e)})^t (0, \dots, 0, \tilde{U}_{(j)}^{(e)})^t.
\end{aligned} 
\quad (36) 
$$

**Step 2.** The boundary layer solution is a pair $(\mathbf{U}_j^{[BLO_l]}, p_j^{[BLO_l]})$ that solves the following Stokes system in $\Omega_l, l = 1, \dots, N_1$:

$$ 
\begin{aligned}
& -\nu \Delta_y \mathbf{U}_j^{[BLO_l]} + \nabla_y P_j^{[BLO_l]} = \mathbf{f}_j^{[REGO_l]} + \mathbf{f}_j^{[BLO_l]}, \\ 
& \mathrm{div}_y \mathbf{U}_j^{[BLO_l]} = h_j^{[REGO_l]}, \\ 
& \mathbf{U}_j^{[BLO_l]}|_{\partial \Omega_l} = 0, \quad j=0, \dots, J,
\end{aligned} 
\quad (37) 
$$

where

$$ 
\begin{aligned}
\mathbf{f}_j^{[REGO_l]}(y^{(e)}) = & -\sum_{e:O_l \in e} \left\{ -\nu \Delta_y \left[ \zeta \left( \frac{y_n^{(e)}}{3r} \right) \mathbf{U}_j^{[e]}(y^{(e)}) \right] \right. \\ 
& \left. + \nabla_y \left[ \zeta \left( \frac{y_n^{(e)}}{3r} \right) (-s_j^{(e)} y_n^{(e)}) \right] \right. \\ 
& \left. + \sum_{p+r=j-1} \zeta \left( \frac{y_n^{(e)}}{3r} \right) (\mathbf{U}_p^{[e]}(y^{(e)}) \cdot \nabla_y) \left( \zeta \left( \frac{y_n^{(e)}}{3r} \right) \mathbf{U}_r^{[e]}(y^{(e)}) \right) \right. \\ 
& \left. + \hat{a}_{j+1}^{(e)} \nabla_y \zeta \left( \frac{y_n^{(e)}}{3r} \right) \right\}
\end{aligned} 
\quad (38) 
$$
(for $j = J$, the coefficient $\hat{a}_{j+1}^{(e)}(t)$ is omitted),

$$ 
\begin{aligned}
\mathbf{f}_j^{[BLO_l]}(y^{(e)}) = & - \sum_{e:O_l \in e} \left\{ \sum_{p+r=j-1} \zeta \left( \frac{y_n^{(e)}}{3r} \right) (\mathbf{U}_p^{[e]}(y^{(e)\prime}) \cdot \nabla_y) \mathbf{U}_r^{[BLO_l]}(y) \right. \\
& \left. + \sum_{p+r=j-1} (\mathbf{U}_p^{[BLO_l]}(y) \cdot \nabla_y) \left( \zeta \left( \frac{y_n^{(e)}}{3r} \right) \mathbf{U}_r^{[e]}(y^{(e)\prime}) \right) \right\} \\
& - \sum_{p+r=j-1} (\mathbf{U}_p^{[BLO_l]}(y) \cdot \nabla_y) \mathbf{U}_r^{[BLO_l]}(y),
\end{aligned} 
\quad (39) $$

$$ h_j^{[REGO_l]}(y,t) = - \sum_{e:O_l \in e} \mathrm{div}_y \left( \zeta \left( \frac{y_n^{(e)}}{3r} \right) \mathbf{U}_j^{[e]}(y^{(e)\prime},t) \right). \quad (40) $$

Here the sum $\sum_{e:O_l \in e}$ is taken over all the edges $e$ that have ends in the node $O_l$, and the terms of the sum are extended by zero out of cylinders $\Pi_e^{(e)}$ and by convention, while the terms with the negative subscripts $j$ are equal to zero.

First, we find the couple $(\mathbf{U}_j^{[BLO_l]}, \hat{P}_j^{[BLO_l]})$, which is the solution to the same Equation (37) without the last term in the definition of $\mathbf{f}_j^{[REGO_l]}$ (see Equation (38)). Using the results of Theorems 7.1 and 7.2 in [5], it can be proven by induction that $\mathbf{U}_j^{[BLO_l]}$ exponentially tends to zero as $|y| \to +\infty$, while the corresponding pressure function $\hat{P}_j^{[BLO_l]}$ stabilizes in outlets to infinity to some constants $\hat{a}_j^{[BLO_l,e]}$ in the sense of Equation (34); 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 to zero for the outlet corresponding to the selected edge $e_s$.

Then we solve the problems in half-cylinders $\Omega_l, l = N_1 + 1, \dots, N$:

$$ 
\begin{aligned}
& -\nu \Delta_y \mathbf{U}_j^{[BLO_l]} + \nabla_y \hat{P}_j^{[BLO_l]} = \mathbf{f}_j^{[REGO_l]} + \mathbf{f}_j^{[BLO_l]}, \\
& \mathrm{div}_y \mathbf{U}_j^{[BLO_l]} = 0, \\
& \mathbf{U}_j^{[BLO_l]}|_{\partial \Omega_l \setminus \{y_n=0\}} = 0, \\
& \mathbf{U}_j^{[BLO_l]}|_{y_n=0} = 0,
\end{aligned} 
\quad (41) $$

$$ 
\begin{aligned}
& \left( -\nu \frac{\partial \mathbf{U}_j^{[BLO_l]} \cdot \mathbf{n}}{\partial \mathbf{n}} + \hat{P}_j^{[BLO_l]} \right) \Big|_{y_n=0} \\
& = -\frac{1}{2} \sum_{p+r=j-1} (\mathbf{U}_p^{[e]}(y^{(e)\prime}) + \mathbf{U}_p^{[BLO_l]}|_{y_n=0}) (\mathbf{U}_r^{[e]}(y^{(e)\prime}) \\
& \qquad + \mathbf{U}_r^{[BLO_l]}|_{y_n=0}), \quad j = 0, \dots, J,
\end{aligned} 
\quad (42) $$

where $\frac{\partial \mathbf{U}_j^{[BLO_l]} \cdot \mathbf{n}}{\partial \mathbf{n}}$ is fictive.

$$ \mathbf{f}_j^{[REGO_l]}(y^{(e)}) = 0, \quad (42) $$

$$ 
\begin{aligned}
\mathbf{f}_j^{[BLO_l]}(y^{(e)}) = & - \sum_{p+r=j-1} (\mathbf{U}_p^{[e]}(y^{(e)\prime}) \cdot \nabla_y) \mathbf{U}_r^{[BLO_l]}(y) \\
& - \sum_{p+r=j-1} (\mathbf{U}_p^{[BLO_l]}(y) \cdot \nabla_y) \mathbf{U}_r^{[e]}(y^{(e)\prime})
\end{aligned} 
\quad (43) $$

$$ - \sum_{p+r=j-1} (\mathbf{U}_p^{[BLO_l]}(y) \cdot \nabla_y) \mathbf{U}_r^{[BLO_l]}(y). $$

The pressure here $\hat{P}_j^{[BLO_l]}$ tends to a constant $\hat{a}_j^{[BLO_l,e]}$.

If $j = J$, then the right-hand side of the boundary condition is replaced by
$$-\frac{1}{2} \sum_{J-1 \le p+r \le 2J} (\mathbf{U}_p^{[e]}(y^{(e)\prime}) + \mathbf{U}_p^{[BLO_l]}|_{y_n=0})(\mathbf{U}_r^{[e]}(y^{(e)\prime}) + \mathbf{U}_r^{[BLO_l]}|_{y_n=0}), \quad j = 0, \dots, J. \quad (44)$$

**Step 3.** Solve the conductivity problem on the graph for the function $p_{j+1}^{(e)}$ ($j < J$):

$$-\kappa_e \frac{\partial^2 p_{j+1}^{(e)}}{\partial x_n^{(e)2}}(x_n^{(e)}) = 0, \quad x_n^{(e)} \in (0, |e|), \\ -\sum_{e:O_l \in e, e \neq e_s} \kappa_e \frac{\partial p_{j+1}^{(e)}}{\partial x_n^{(e)}}(0) = 0, \quad l = 1, \dots, N_1, \\ p_{j+1}^{(e)}(0) = \hat{a}_j^{[BLO_l,e]}, \quad l = N_1+1, \dots, N, \\ p_{j+1}^{(e)}(0) - p_{j+1}^{(e_s)}(0) = \hat{a}_j^{[BLO_l,e]}, \quad \forall e \subset B_l, e \neq e_s,$$

where the local coordinates $x^{(e)}$ are redefined so that all of them have the same origin $O_l$.

**Step 4.** Finally, we find the pressure $P_j^{[BLO_l]}(y)$ in the boundary layer problem (Equations (37) and (38)) for $l = 1, \dots, N_1$:

$$P_j^{[BLO_l]}(y) = \hat{P}_j^{[BLO_l]}(y) - \sum_{e:O_l \in e, e \neq e_s} \zeta\left(\frac{y_n^{(e_{\alpha m})}}{3r}\right) \hat{a}_j^{[BLO_l,e]},$$

and $P_j^{[BLO_l]}(y)$ in the boundary layer problem (Equations (41) and (42)) for $l = N_1 + 1, \dots, N$:

$$P_j^{[BLO_l]}(y) = \hat{P}_j^{[BLO_l]}(y) - \hat{a}_j^{[BLO_l,e]}.$$

This step finalizes the passage from $j$ to $j+1$.

## 5. Residual

Consider the asymptotic expansion $(\mathbf{u}^{(J)}, p^{(J)})$ of order $J$ (see Equations (22) and (24)). By the construction

$$\mathbf{u}^{(J)} \in W^{2,2}(B_{\varepsilon}), \nabla p^{(J)} \in L^2(B_{\varepsilon}). \quad (45)$$

Moreover, $\|\mathbf{u}^{(J)}\|_{L^4(B_{\varepsilon})}^4 \le c\varepsilon^{(n-1)/4}$. Indeed, the Poiseuille part of $\mathbf{u}^{(J)}$ satisfies this estimate. The $\|\cdot\|_{L^4(B_{\varepsilon}^{(i)})}$-norm of the boundary layer functions in each bundle $B_{\varepsilon}^{(i)}$ can be estimated by the $L^4$-norm in the unbounded dilated domain $\Omega_i$ multiplied by $\varepsilon^{n/4}$. Taking into consideration an exponential decay of the boundary layers, we obtain the desired estimate.

Put $\mathcal{L}(\mathbf{u}, p) = -\nu\Delta\mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{u} + \nabla p$. Let us calculate $\mathcal{L}(\mathbf{u}^{(J)}, p^{(J)})$ in a half-bundle $HB_{O_l}, l = 1, \dots, N_1$. We obtain

$$\mathbf{f}^{(J)}(x) = \mathcal{L}(\mathbf{u}^{(J)}, p^{(J)}) \\ = \sum_{J+1 \le j \le 2J} \varepsilon^{j-2} \left( \sum_{e:O_l \in e} \sum_{\alpha+\beta=j-1} \left( \mathbf{U}_{(\alpha)}^{(e)} \zeta\left(\frac{y_n^{(e)}}{3r}\right) \cdot \nabla_y \right) \mathbf{U}_{(\beta)}^{(e)} \zeta\left(\frac{y_n^{(e)}}{3r}\right) \right. \\ + \sum_{e:O_l \in e} \left[ \sum_{p+r=j-1} \zeta\left(\frac{y_n^{(e)}}{3r}\right) \left( \mathbf{U}_{(p)}^{(e)}(y^{(e)\prime}) \cdot \nabla_y \right) \mathbf{U}_r^{[BLO_l]}(y) \right. \\ \left. + \sum_{p+r=j-1} \left( \mathbf{U}_p^{[BLO_l]}(y,t) \cdot \nabla_y \right) \left( \zeta\left(\frac{y_n^{(e)}}{3r}\right) \mathbf{U}_r^{(e)}(y^{(e)\prime}) \right) \right)$$
$$ + \sum_{p+r=j-1} (\mathbf{U}_p^{[BLO_l]}(y,t) \cdot \nabla_y) \mathbf{U}_r^{[BLO_l]}(y) ] + \varepsilon^{J-2} \sum_{e:O_l \in e} \hat{a}_j^{[BLO_{l,e}]} \nabla_y \zeta \left( \frac{y_n^{(e)}}{3r} \right) - \left\{ \mathcal{L} \left( \zeta \left( \frac{|x-O_l|}{|e|_{min}} \right) \mathbf{U}^{[BLO_{l,J}]}(y,t), \zeta \left( \frac{|x-O_l|}{|e|_{min}} \right) P^{[BLO_{l,J}]}(y,t) \right) \chi(x) \right\}. $$

Here $y = \frac{x - O_l}{\varepsilon}$; $y^{(e)} = \frac{x^{(e)}}{\varepsilon}$; $\chi = \chi_{\text{supp}}(1 - \zeta(\frac{|x-O_l|}{|e|_{min}}))$ is the characteristic function of the set $\text{supp}(1 - \zeta(\frac{|x-O_l|}{|e|_{min}}))$. As before, the terms of the sums $\sum_{e:O_l \in e}$ are extended by zero out of cylinders $\Pi_{\varepsilon}^{(e)}$.

Here the first four lines come from the inertial term, and they contain all the combinations of $\mathbf{U}_{\beta}^{(e)}$ and $\mathbf{U}_{\beta}^{[BLO_l]}$ that have an order higher than $J-2$, and the next line 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 $\zeta(\frac{|x-O_l|}{|e|_{min}})$. Notice that terms appearing in this last line exponentially vanish because in the set $\text{supp}(1 - \zeta(\frac{|x-O_l|}{|e|_{min}}))$ (where $\chi \neq 0$), the order of this term in $L^2$-norm is $O(e^{-c_1/\varepsilon})$ with some positive constant $c_1$ (see the Appendix in [5]). From the obtained formulas, it follows that

$$ \|\mathbf{f}^{(J)}\|_{L^2(B_{\varepsilon})} = \|\mathcal{L}(\mathbf{u}^{(J)}, p^{(J)})\|_{L^2(B_{\varepsilon})} = O(\varepsilon^{J-2}). \quad (46) $$

In the vertex-associated cylinders $B_{O_l}$, $l = N_1 + 1, \dots, N$, the residual is simpler: it is without the factor $\zeta(\frac{y_n^{(e)}}{3r})$.

Let us calculate the divergence of $\mathbf{u}^{(J)}$. In any half-bundle, we have

$$ \text{div}\mathbf{u}^{(J)} = -\nabla\zeta\left(\frac{|x - O_l|}{|e|_{min}}\right) \cdot \mathbf{U}^{[BLO_l, J]}(y) = h^{(J)}(y). $$

Obviously, $h^{(J)} \in W^{1,2}(B_{\varepsilon})$. Since the support of the function $\nabla\zeta(\frac{|x-O_l|}{|e|_{min}})$ belongs to the middle third of every cylinder, the relations there

$$ \|h^{(J)}\|_{W^{1,2}(B_{\varepsilon})} = O(e^{-c_2/\varepsilon}) \quad (47) $$

hold for some $c_2 > 0$.

Finally, the boundary conditions are satisfied with the residual $\varepsilon^{J-1}\hat{a}_j^{[BLO_{l,e}]}$ on $\gamma_{\varepsilon}^l$. This residual appears as a result of the subtraction of the constant $\hat{a}_j^{[BLO_{l,e}]}$ from the boundary layer pressure $P_j^{[BLO_l]}(y)$ in Step 4 of the algorithm. For all $j < J$, it is compensated by the gaps of the pressure in the problem on the graph, but for $j = J$, it remains as a residual.

It is easy to see that

$$ \int_{B_{\varepsilon}} h^{(J)}(y) dy = 0. $$

Therefore, by Lemma 3.7 in [17], there exists a vector field $\mathbf{w}^{(J)} \in \hat{W}^{1,2}(B_{\varepsilon})$ with such that $\text{div}\mathbf{w}^{(J)} = -h^{(J)}$. Moreover, the estimates

$$ \|\mathbf{w}^{(J)}\|_{W^{1,2}(B_{\varepsilon})} \leq \varepsilon^{-1}c\|h^{(J)}\|_{L^2(B_{\varepsilon})} \quad (48) $$

hold.
Set $\tilde{\mathbf{u}}^{(J)} = \mathbf{u}^{(J)} + \mathbf{w}^{(J)}$. Then $\mathrm{div}\tilde{\mathbf{u}}^{(J)} = 0$, $\tilde{\mathbf{u}}^{(J)}$ satisfies the boundary conditions with the residual $-\varepsilon^{J-1}\hat{a}_j^{[BLO_l, e]}$ on $\gamma_\varepsilon^l$, and because of Equation (47), we have

$$ \|\mathbf{f}_1^{(J)}\|_{L^2(B_\varepsilon)} = O(\varepsilon^{J-2}), \qquad (49) $$

where $\mathbf{f}_1^{(J)} = \mathcal{L}(\tilde{\mathbf{u}}^{(J)}, p^{(J)})$.

## 6. Error Estimate

**Theorem 3.** The following error estimate

$$ \|\mathbf{u} - \tilde{\mathbf{u}}^{(J)}\|_{W^{1,2}(B_\varepsilon)} = O(\varepsilon^{J+(n-1)/2}) \qquad (50) $$

holds.

**Proof.** Let $\mathbf{v} = \mathbf{u} - \tilde{\mathbf{u}}^{(J)}$. Then the integral identity

$$ \begin{aligned} & \nu \int_{B_\varepsilon} \nabla \mathbf{v} : \nabla \eta \, dx + \int_{B_\varepsilon} (\mathbf{v} \cdot \nabla) \mathbf{v} \cdot \eta \, dx + \int_{B_\varepsilon} (\mathbf{v} \cdot \nabla) \tilde{\mathbf{u}}^{(J)} \cdot \eta \, dx \\ & \qquad + \int_{B_\varepsilon} (\tilde{\mathbf{u}}^{(J)} \cdot \nabla) \mathbf{v} \cdot \eta \, dx - \int_{B_\varepsilon} (\eta \cdot \nabla) \mathbf{v} \cdot \mathbf{v} \, dx \\ & \qquad - \int_{B_\varepsilon} (\eta \cdot \nabla) \mathbf{v} \cdot \tilde{\mathbf{u}}^{(J)} \, dx - \int_{B_\varepsilon} (\eta \cdot \nabla) \tilde{\mathbf{u}}^{(J)} \cdot \mathbf{v} \, dx \\ & = \varepsilon^{J-1} \sum_{l=N_1+1}^{N} \hat{a}_j^{[BLO_l, e]} \int_{\gamma_\varepsilon^l} \eta_n \, dx' - \int_{B_\varepsilon} \mathbf{f}_1^{(J)} \cdot \eta \, dx \end{aligned} $$

holds for every $\eta \in \hat{\gamma}_\varepsilon^{1,2}(B_\varepsilon)$.

Taking $\eta = \mathbf{v}$ and integrating by parts, we obtain

$$ \begin{aligned} & \nu \int_{B_\varepsilon} |\nabla \mathbf{v}|^2 \, dx = \varepsilon^{J-1} \sum_{l=N_1+1}^{N} \hat{a}_j^{[BLO_l, e]} \int_{\gamma_\varepsilon^l} \mathbf{v}_n \, dx' \\ & + \int_{B_\varepsilon} (\mathbf{v} \cdot \nabla) \mathbf{v} \cdot \tilde{\mathbf{u}}^{(J)} \, dx - \int_{B_\varepsilon} (\tilde{\mathbf{u}}^{(J)} \cdot \nabla) \mathbf{v} \cdot \mathbf{v} \, dx - \int_{B_\varepsilon} \mathbf{f}_1^{(J)} \cdot \mathbf{v} \, dx. \end{aligned} \qquad (51) $$

From Lemma 3, it follows that

$$ \left| \varepsilon^{J-1} \sum_{l=N_1+1}^{N} \hat{a}_j^{[BLO_l, e]} \int_{\gamma_\varepsilon^l} \mathbf{v}_n \, dx' \right| \le c\varepsilon^{J-1+n/2} \sum_{l=N_1+1}^{N} \hat{a}_j^{[BLO_l, e]} \|\nabla \mathbf{v}\|_{L^2(B_\varepsilon)} \cdot $$

Using the Hölder inequality, Equations (2) and (3), and the estimate $\|\tilde{\mathbf{u}}^{(J)}\|_{L^4(B_\varepsilon)} \le c\varepsilon^{(n-1)/4}$, we obtain

$$ \begin{aligned} & \left| \int_{B_\varepsilon} (\mathbf{v} \cdot \nabla) \mathbf{v} \cdot \tilde{\mathbf{u}}^{(J)} \, dx \right|, \left| \int_{B_\varepsilon} (\tilde{\mathbf{u}}^{(J)} \cdot \nabla) \mathbf{v} \cdot \tilde{\mathbf{u}}^{(J)} \, dx \right| \le \|\mathbf{v}\|_{L^4(B_\varepsilon)} \|\nabla \mathbf{v}\|_{L^2(B_\varepsilon)} \|\tilde{\mathbf{u}}^{(J)}\|_{L^4(B_\varepsilon)} \\ & \le c\varepsilon^{\alpha+\frac{n-1}{4}} \|\nabla \mathbf{v}\|_{L^2(B_\varepsilon)}^2 \le c\varepsilon^{\frac{3}{4}} \|\nabla \mathbf{v}\|_{L^2(B_\varepsilon)}^2, \end{aligned} $$

where $\alpha = 1/2$ for $n=2$ and $\alpha = 1/4$ for $n=3$. Moreover,

$$ \left| \int_{B_\varepsilon} \mathbf{f}_1^{(J)} \cdot \mathbf{v} \, dx \right| \le \varepsilon C_{PF} \|\mathbf{f}_1^{(J)}\|_{L^2(B_\varepsilon)} \|\nabla \mathbf{v}\|_{L^2(B_\varepsilon)} $$

where $\varepsilon C_{PF}$ is the Poincaré–Friedrichs constant for the domain $B_\varepsilon$.
From these estimates and Equation (51), we obtain

$$ \nu \int_{B_{\epsilon}} |\nabla \mathbf{v}|^2 dx \le c\epsilon^{\frac{2(J-1)+n}{2}} \sum_{l=N_1+1}^{N} \hat{a}_J^{[BLO_l, \epsilon]} \|\nabla \mathbf{v}\|_{L^2(B_\epsilon)} \\ +c\epsilon^{\frac{3}{4}} \|\nabla \mathbf{v}\|_{L^2(B_\epsilon)}^2 + \epsilon C_{PF} \|\mathbf{f}_1^{(J)}\|_{L^2(B_\epsilon)} \|\nabla \mathbf{v}\|_{L^2(B_\epsilon)}. $$

Hence

$$ \|\nabla \mathbf{v}\|_{L^2(B_\epsilon)} \le \frac{C}{\nu - c\epsilon^{3/4}} \left( \epsilon^{\frac{2(J-1)+n}{2}} \sum_{l=N_1+1}^{N} \hat{a}_J^{[BLO_l, \epsilon]} + \epsilon \|\mathbf{f}_1^{(J)}\|_{L^2(B_\epsilon)} \right) \\ \le \frac{C}{\nu - c\epsilon^{3/4}} \left( \epsilon^{\frac{2(J-1)+n}{2}} \sum_{l=N_1+1}^{N} \hat{a}_J^{[BLO_l, \epsilon]} + \epsilon^{J-1} \right). $$

If $\nu - c\epsilon^{3/4} > \nu/2$, then

$$ \|\mathbf{v}\|_{W^{1,2}(B_\epsilon)} = \|\mathbf{u} - \tilde{\mathbf{u}}^{(J)}\|_{W^{1,2}(B_\epsilon)} = O(\epsilon^{J-1}). \qquad (52) $$

Now, by evaluating the norm of the difference $\mathbf{u}^{(J)}$ and $\mathbf{u}^{(J+2)}$, we obtain:

$$ \|\tilde{\mathbf{u}}^{(J)} - \tilde{\mathbf{u}}^{(J+2)}\|_{W^{1,2}(B_\epsilon)} = O(\epsilon^{J+(n-1)/2}). $$

Replacing $J$ with $J+2$ in Equation (52), we obtain:

$$ \|\mathbf{u} - \tilde{\mathbf{u}}^{(J+2)}\|_{W^{1,2}(B_\epsilon)} = O(\epsilon^{J+1}). $$

So, the triangle inequality gives Equation (50). □

## 7. Conclusions

The main result of the paper is the construction of the asymptotic expansion of a weak solution of the stationary Navier-Stokes equations in a thin tube structure with the Bernoulli boundary conditions for the inflows and outflows. The existence and uniqueness of the solution is proved. The dependence of the stability estimate on the small parameter is addressed. We proved the error estimate. It allows evaluating the limitations of the theoretical predictions of the asymptotic theory and introducing a new numerical strategy for computations of the Navier-Stokes equations in thin tube structures.

**Author Contributions:** Conceptualization, R.J., G.P. and K.P.; methodology, R.J., G.P. and K.P.; formal analysis, R.J., G.P. and K.P.; investigation, R.J., G.P. and K.P.; resources, R.J., G.P. and K.P.; writing—original draft preparation, R.J., G.P. and K.P.; writing—review and editing, R.J., G.P. and K.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research received an external funding from the LABEX MILYON, Russian Science Foundation and European social Fund (see Acknowledgement).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** R. Juodagalvytė was supported by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency. G. Panasenko was supported by the Russian Science Foundation grant 19-11-00033 operated by the Moscow Power Engineering Institute and K. Pileckas by the European Social Fund (project No 09.3.3-LMT-K-712-17-0003) under grant agreement with the Research Council of Lithuania (LMTLT).

**Conflicts of Interest:** The authors declare no conflict of interest.
## Appendix A. Stokes Equation in a Half-Cylinder with Neumann's Condition on the Base and No-Slip Condition on the Lateral Boundary

Let $\Omega$ be a half-cylinder $\omega \times (0, +\infty)$, where $\omega$ is a bounded domain in $\mathbb{R}^{n-1}$ with a Lipschitz boundary. $\Gamma$ denotes the lateral boundary $\partial\omega \times (0, +\infty)$, and $\gamma$ denotes the base $\omega \times \{0\}$. Consider the stationary Stokes problem

$$ egin{cases} 
-v\Delta\mathbf{u} + \nabla p = \mathbf{f} + \sum_{m=1}^{n} \frac{\partial \mathbf{f}_m}{\partial x_m}, \\ 
\operatorname{div}\mathbf{u} = 0, \\ 
\mathbf{u}(x)|_{\Gamma} = 0, \\ 
p(x)|_{\gamma} = \psi(x'), \\ 
\mathbf{u}_{\tau}(x)|_{\gamma} = 0. 
\end{cases} 
\quad (A1) $$

Define $J_{\Gamma,0} = \{\eta \in W^{1,2}(\Omega) : \operatorname{div}\eta = 0, \eta|_{\Gamma} = 0, \eta_{\tau}|_{\gamma} = 0\}$. Assume that $\mathbf{f}, \mathbf{f}_m \in L^2(\Omega)$ and $\psi \in L^2(\gamma)$. By a weak solution of problem (A1), we understand a vector field $\mathbf{u} \in J_{\Gamma,0}$ satisfying the integral identity

$$ v \int_{\Omega} \nabla \mathbf{u} : \nabla \eta \,dx = \int_{\gamma} \psi(x') \eta_n(x') \,dx' + \int_{\Omega} \mathbf{f} \cdot \eta \,dx - \sum_{m=1}^{n} \int_{\Omega} \mathbf{f}_m \cdot \frac{\partial \eta}{\partial x_m} \,dx \quad (A2) $$

for every vector field $\eta \in J_{\Gamma,0}$.

**Theorem A1.** Assume that $\mathbf{f}, \mathbf{f}_m \in L^2(\Omega)$ and $\psi \in L^2(\gamma)$. Then there exists a unique weak solution $\mathbf{u}$ of problem (A1). It satisfies the estimate

$$ \|\mathbf{u}\|_{W^{1,2}(\Omega)} \le C \left( \|\mathbf{f}\|_{L^2(\Omega)} + \sum_{m=1}^{n} \|\mathbf{f}_m\|_{L^2(\Omega)} + \|\psi\|_{L^2(\gamma)} \right) \quad (A3) $$

with a constant $C$ independent of $\mathbf{f}, \mathbf{f}_m$ and $\psi$.

**Proof.** The proof of this theorem is standard: the right-hand side of Equation (A2) is considered as a linear continuous functional on the space of functions $\eta$ belonging to $J_{\Gamma,0}$. This space is supplied with the inner product $[\mathbf{u}, \eta] = v \int_{\Omega} \nabla \mathbf{u} \cdot \nabla \eta \,dx$. Using the trace theorem and the Poincaré-Friedrichs inequality for functions from $J_{\Gamma,0}$, we apply the Riesz representation theorem and obtain the existence and uniqueness of the solution. Then taking $\eta = \mathbf{u}$ and using again the trace theorem and Poincaré-Friedrichs inequality, we obtain estimate (A3). $\square$

Let us define in $\Omega$ weighted function spaces. Denote

$$ E_{\beta}(x) = \exp(2\beta x_n). \quad (A4) $$

Denote by $W_{\beta}^{l,2}(\Omega)$, $l \ge 0$ the space of functions obtained as the closure of $C_0^{\infty}(\Omega)$ in the norm

$$ \|u\|_{W_{\beta}^{l,2}(\Omega)} = \left( \sum_{|\alpha|=0}^{l} \int_{\Omega} E_{\beta}(x) |D^{\alpha}u(x)|^2 dx \right)^{1/2} $$

and set $W_{\beta}^{0,2}(\Omega) = \mathcal{L}_{\beta}^2(\Omega)$. Notice that for $\beta > 0$, elements of the space $W_{\beta}^{l,2}(\Omega)$ exponentially vanish as $x_n \to \infty$.

Denote $\Omega_{\delta} = \{x \in \Omega : x_n > \delta\}$.

There holds the following theorem.

**Theorem A2.** Assume that $\mathbf{f}, \mathbf{f}_m \in \mathcal{L}_{\beta}^2(\Omega)$, $\beta > 0$. If $\beta$ is sufficiently small, then the weak solution $\mathbf{u}$ of problem (A1) belongs to the space $W_{\beta}^{1,2}(\Omega)$.
Moreover, if $\partial\omega \in C^2$ and $f_m = 0$, then for any $\delta > 0$, $\mathbf{u} \in W_{\beta}^{2,2}(\Omega_{\delta})$, and there exists a function $p \in L_{\text{loc}}^2(\Omega)$ with $\nabla p \in \mathcal{L}_{\beta}^2(\Omega_{\delta})$ such that the pair $(\mathbf{u}, p)$ satisfies Equation (A1) almost everywhere in $\Omega_{\delta}$. There exists a constant $\hat{a}$ such that $\lim_{x \in \Omega, |x| \to \infty} p(x) = \hat{a}$ in the sense

$$ \int_{\Omega_{\delta}} \exp\{2\beta_1 x_n\} |p(x) - \hat{a}|^2 dx < \infty \quad \forall \beta_1 \in (0, \beta). \quad (A5) $$

This assertion is a corollary of Theorems A.1 and A.2 and Proposition A.1 of [5], see also [18,19]. The regularity of the solution in $\Omega_{\delta}$ needed for the proof follows from ADN estimates (see [20]).

## References

1. Panasenko, G. Asymptotic expansion of the solution of Navier-Stokes equation in a tube structure. *Comptes Rendus Acad. Sci. Ser. IIB Mech. Phys. Astron.* **1998**, 326, 867–872. [CrossRef]

2. Galdi, P.; Rannacher, R.; Robertson, A.M.; Turek, S. *Hemodynamical Flows, Modeling, Analysis and Simulation*; Oberwolfach Seminars; Birkhäuser: Basel, Switzerland, 2008; Volume 37.

3. Formaggia, L.; Lamponi, D.; Quarteroni, A. One dimensional models for blood flow in arteries. *J. Eng. Math.* **2003**, 47, 251–276. [CrossRef]

4. Panasenko, G. Partial asymptotic decomposition of domain: Navier-Stokes equation in tube structure. *Comptes Rendus Acad. Sci. Ser. IIB Mech. Phys. Astron.* **1998**, 326, 893–898. [CrossRef]

5. Panasenko, G.; Pileckas, K. Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure. I. The case without boundary layer-in-time. *Nonlinear Anal. Theory Methods Appl.* **2015**, 122, 125–168. [CrossRef]

6. Conca, C.; Murat, F.; Pironneau, O. The Stokes and Navier-Stokes equations with boundary conditions involving pressure. *Jpn. J. Math.* **1994**, 20, 2. [CrossRef]

7. Marušić-Paloka, E. Rigorous justification of the Kirchhoff law for junction of thin pipes with viscous fluid. *Asymptot. Anal.* **2003**, 33, 51–66.

8. Kazhikhov, A.V. Solvability of some unilateral boundary value problems for Navier-Stokes equations. In *Dynamics of Continuous Media*; Institute of Hydrodynamics of Siberian Department of Academy of Sciences of USSR: Novosibirsk, Russia, 1974; Volume 16, pp. 5–34.

9. Korobkov, M.V.; Pileckas, K.; Russo, R. Solvability in a finite pipe of steady-state Navier-Stokes equations with boundary conditions involving Bernoulli pressure. *Calc. Var. Partial. Differ. Equ.* **2020**, 59, 32. [CrossRef]

10. Baranovskii, E.S. A novel 3D model for non-Newtonian fluid flows in a pipe network. *Math. Methods Appl. Sci.* **2021**, 44, 3827–3839. [CrossRef]

11. Panasenko, G. *Multi-Scale Modeling for Structures and Composites*; Springer: Dordrecht, The Netherlands, 2005.

12. Juodagalvytė, R.; Panasenko, G.; Pileckas, K. Time periodic Navier-Stokes equations in a thin tube structure. *Bound. Value Probl.* **2020**, 2020, 28. [CrossRef]

13. Egloffe, A.-C. Study of Some Inverse Problems for the Stokes System. Application to the Lungs. Ph.D. Thesis, l'Université Pierre et Marie Curie, Paris, France, 2012.

14. Ladyzhenskaya, O.A. *The Mathematical Theory of Viscous Incompressible Fluid*; Gordon and Breach: New York, NY, USA, 1969.

15. Heywood, J.G.; Rannacher, R.; Turek, S. Artificial boundaries and flux and pressure conditions for the incompressible Navier-Stokes equations. *Int. J. Numer. Methods Fluids* **1996**, 22, 325–352. [CrossRef]

16. Panasenko, G.; Pileckas, K. Flows in a tube structure: Equation on the graph. *J. Math. Phys.* **2014**, 55, 081505. [CrossRef]

17. Panasenko, G.; Pileckas, K. Divergence equation in thin-tube structure. *Appl. Anal.* **2015**, 94, 1450–1459. [CrossRef]

18. Pileckas, K. On the nonstationary linearized Navier-Stokes problem in domains with cylindrical outlets to infinity. *Math. Ann.* **2005**, 332, 395–419. [CrossRef]

19. Pileckas, K. Navier-Stokes system in domains with cylindrical outlets to infinity. Leray's problem. In *Handbook of Mathematical Fluid Dynamics*; Elsevier: Amsterdam, The Netherlands, 2007; Chapter 8, Volume 4, pp. 445–647.

20. Agmon, S.; Douglis, A.; Nirenberg, L. Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions. *I. Commun. Pure Appl. Math.* **1959**, 12, 623–727. [CrossRef]