# Steady state non-Newtonian flow with strain rate dependent viscosity in thin tube structure with no slip boundary condition

Grigory Panasenko, Konstantin Pileckas, Bogdan Vernescu

► To cite this version:

Grigory Panasenko, Konstantin Pileckas, Bogdan Vernescu. Steady state non-Newtonian flow with strain rate dependent viscosity in thin tube structure with no slip boundary condition. Mathematical Modelling of Natural Phenomena, 2022, 17, pp.18. 10.1051/mmnp/2022005. hal-04134122

HAL Id: hal-04134122

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

Submitted on 13 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
# STEADY STATE NON-NEWTONIAN FLOW WITH STRAIN RATE DEPENDENT VISCOSITY IN THIN TUBE STRUCTURE WITH NO SLIP BOUNDARY CONDITION

GRIGORY PANASENKO<fnref n="1" />,*, KONSTANTIN PILECKAS<fnref n="2" />
AND BOGDAN VERNESCU<fnref n="3" />

**Abstract.** The steady state non-Newtonian flow, with strain rate dependent viscosity in a thin tube structure, with no slip boundary condition, is considered. Applying the Banach fixed point theorem we prove the existence and uniqueness of a solution. An asymptotic approximation is constructed and justified by an error estimate.

**Mathematics Subject Classification.** 35Q35, 76D07.

Received October 22, 2021. Accepted January 24, 2022.

## 1. INTRODUCTION

The asymptotic behavior of solutions of partial differential equations in thin domains is extensively studied in a vast mathematical literature. In particular, the tube structures, introduced in [20], are considered as a geometrical model of a blood vessels network. Viscous flows were studied in such domains and for the steady state Navier-Stokes equations an asymptotic expansion of the solution was constructed. The non-stationary Navier-Stokes equations in such a domain were studied in [24]. A non-Newtonian flow (Bingham flow) in a network of thin pipes was studied mathematically in [4]. However the asymptotic behavior is described there neglecting the boundary layer functions. The power-law rheology in a thin tube structures was studied in [16], where the first order asymptotic approximation with boundary layers was constructed (see [15] for the existence and uniqueness of a solutions to boundary layer problems). In the present paper we consider the flow with a strain rate dependent viscosity, taking into consideration the boundary layers. We will construct an asymptotic expansion of the solution by an iterative algorithm. The estimates are obtained for higher orders of asymptotic approximations. The main application of these results is modeling for the blood flow in a network of thin vessels.

The leading term of asymptotic expansion of the pressure is described by a one-dimensional elliptic non-linear problem on the graph. The pressure is a linear function on each edge of the graph. At the junctions (nodes) the Kirchhoff like conditions are set. On the other hand, one-dimensional models derived from the conservation laws

<footnote>*Keywords and phrases:* Non-Newtonian flow, strain rate dependent viscosity, asymptotic approximation, quasi-Poiseuille flows, boundary layers, method of asymptotic partial decomposition.</footnote>

<footnote n="1">¹ University of Lyon, UJM, Institute Camille Jordan UMR CNRS 5208, 23 rue P. Michelon, 42023, Saint-Etienne, France and Institute of Applied Mathematics, Vilnius University, Naugarduko Str., 24, Vilnius 03225, Lithuania.</footnote>

<footnote n="2">² Institute of Applied Mathematics, Vilnius University, Naugarduko Str., 24, Vilnius 03225, Lithuania.</footnote>

<footnote n="3">³ Department of Mathematical Sciences, Worcester Polytechnic Institute, Worcester, MA 01609, USA.</footnote>

<footnote>* Corresponding author: grigory.panasenko@univ-st-etienne.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Graph B.</figcaption></figure>

were introduced in [5]. These models differ from the problem on the graph: they are equations of hyperbolic type. Below in the Conclusion we will discuss the limitations of both approaches, i.e. of the problem on the graph and of the hyperbolic system of equations.

Let $G$ be a bounded domain in $\mathbb{R}^n$. By $L^p(G)$ and $W^{m,p}(G)$, $1 < p < \infty$, $m \ge 1$, we denote the usual Lebesgue and Sobolev spaces, respectively. The norms in $L^p(G)$ and $W^{m,p}(G)$ are indicated by $\| \cdot \|_{L^p(G)}$ and $\| \cdot \|_{W^{m,p}(G)}$, respectively. $W^{m-1/2,2}(\partial G)$ is the space of traces on $\partial G$ of functions from $W^{m,2}(G)$.

By $\dot{W}^{m,p}(G)$ we denote the closure of the set $C_0^\infty(G)$ in the norm $\| \cdot \|_{W^{m,p}(G)}$, where $C_0^\infty(G)$ is the set of all infinitely differentiable functions with compact supports in $G$. More information about these spaces can be found in [1]. Vector-valued functions are denoted by bold letters, and the spaces of scalar and vector-valued functions are not distinguished in notation. More specific functions spaces are introduced in places where they are used.

## 1.1. Thin tube structure

Let us recall the definitions of the tube structure and its graph given in [21].

**Definition 1.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 = O_{ij}O_{kj}$, 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 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 $\mathcal{B} = \bigcup_{j=1}^M e_j$ the union of edges and assume that $\mathcal{B}$ is a connected set (see Fig. 1). The union of all edges having the same end point $O_l$ is called the bundle $\mathcal{B}^{(l)}$.

Let $e$ be some edge, $e = 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_1^{(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_i \tilde{x}_1^{(e)}$ is directed over the ray $[O_jO_i)$.

Below in 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 $\sigma^j \subset \mathbb{R}^{n-1}$ containing the origin $O$ and having $C^4$-smooth boundary $\partial\sigma^j$, $j = 1, \dots, M$. For every edge $e_j = e$ and associated $\sigma^j = \sigma^{(e)}$ we denote by $\Pi_{\varepsilon}^{(e)}$ the cylinder

$$ \Pi_{\varepsilon}^{(e)} = \{ x^{(e)} \in \mathbb{R}^n : x_1^{(e)} \in (0, |e|), \frac{x^{(e)}'}{\varepsilon} \in \sigma^{(e)} \}, $$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Thin tube structure $B_\epsilon$.</figcaption></figure>

where $x^{(e)\prime} = (x_2^{(e)}, \dots, x_n^{(e)})$, $|e|$ is the length of the edge $e$ and $\epsilon > 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 $\epsilon$. We will define as well a semi-infinite dilated cylinder $\Pi_\infty^{(e)} = \{x^{(e)} \in \mathbb{R}^n : x_1^{(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$ be bounded independent of $\epsilon$ domains in $\mathbb{R}^n$; introduce the nodal domains $\omega_\epsilon^j = \{x \in \mathbb{R}^n : \frac{x - O_j}{\epsilon} \in \omega^j\}$.

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

**Definition 1.2.** By a tube structure (see Fig. 2) we call the following domain

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

Suppose that it is a connected set and that the boundary $\partial B_\epsilon$ of $B_\epsilon$ is $C^4$-smooth.

Let $r_1$ be the maximal diameter of domains $\omega_i$, $i = 1, \dots, 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_\epsilon^l \cup \left( \bigcup_{j:O_l \in e_j} \Pi_\epsilon^{(e_j)} \right)$, and by $\Omega_l = \omega^l \cup \left( \bigcup_{j:O_l \in e_j} \Pi_\epsilon^{(e_j)} \right)$ a bundle of dilated cylinders. Denote also $\Omega_l^\epsilon = \{x \in \mathbb{R}^n | x/\epsilon \in \Omega_l\}$.

## 1.2. Formulation of the problem

Let $\nu_0, \lambda > 0$ be positive constants. Let $\nu$ be a bounded $C^3$-smooth function $\mathbb{R}^{n(n+1)/2} \to \mathbb{R}$ such that for all $y \in \mathbb{R}^{n(n+1)/2}$,

$$|\nu(y)| \le A, \quad |\nabla \nu(y)| \le A, \quad |\nabla^2(\nu(y))| \le A, \quad |\nabla^3(\nu(y))| \le A. \qquad (1.1)$$

where $A$ is a positive constant independent of $y$.

Consider in the tube structure $B_\epsilon$ the steady state boundary value problem for the non-Newtonian fluid motion equations

$$ \left\{ \begin{array}{l} -\mathrm{div}((\nu_0 + \lambda\nu(\dot{\gamma}(\mathbf{v}))D(\mathbf{v})) + \nabla p) = 0, \quad x \in B_{\epsilon}, \\ \mathrm{div}\mathbf{v} = 0, \quad x \in B_{\epsilon}, \\ \mathbf{v}|_{\partial B_{\epsilon}} = \epsilon\mathbf{g}, \end{array} \right. \qquad (1.2) $$
where $D(\mathbf{v})$ is the strain rate matrix with the elements $d_{ij} = \frac{1}{2}(\frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i})$, $\dot{\gamma}(\mathbf{v}) = (d_{12}, d_{13}, d_{23}, d_{11}, d_{22}, d_{33})$ if $n = 3$ and $\dot{\gamma}(\mathbf{v}) = (d_{12}, d_{11}, d_{22})$ if $n = 2$.<fnref n="1" />

Assume that the fluid velocity $\mathbf{g}$ at the boundary $\partial B_{\varepsilon}$ has the following structure: $\mathbf{g} = 0$ everywhere on $\partial B_{\varepsilon}$ except for the set $\gamma_{\varepsilon}^{N_1+1}, \dots, \gamma_{\varepsilon}^N$, where $\gamma_{\varepsilon}^j = \partial B_{\varepsilon} \cap \partial \omega_{\varepsilon}^j$, $j = N_1 + 1, \dots, N$, i.e.,

$$ \mathbf{g}(x)|_{\gamma_{\varepsilon}^j} = \mathbf{g}^j \left( \frac{x - O_j}{\varepsilon} \right) \Big|_{\gamma_{\varepsilon}^j}, \quad j = N_1 + 1, \dots, N, \\ \mathbf{g}(x, t)|_{\partial B_{\varepsilon} \setminus (\bigcup_{j=N_1+1}^N \gamma_{\varepsilon}^j)} = 0, $$

where $\mathbf{g}^j \in W^{5/2,2}(\gamma^j)$, $\gamma^j = \varepsilon^{-1}(\gamma_{\varepsilon}^j - O_j)$ are the corresponding dilated parts of the boundary, $\mathbf{g} \in W^{5/2,2}(\partial B_{\varepsilon})$. Assume that

$$ \sum_{j=N_1+1}^{N} \int_{\gamma_{\varepsilon}^j} \mathbf{g}^j \left( \frac{x - O_j}{\varepsilon} \right) \cdot \mathbf{n}(x) dS = 0, \quad (1.3) $$

where $\mathbf{n}$ is the unit outward (with respect to $B_{\varepsilon}$) normal vector to $\gamma_{\varepsilon}^j$.

## 1.3. Main results

The first main result of the paper is the theorem on the existence and uniqueness of the solution: There exists $\lambda_0$ such that for all $\lambda \in (0, \lambda_0)$ problem (1.2) admits a unique weak solution $(\mathbf{u}, p)$ with $\mathbf{u} \in W^{3,2}(B_{\varepsilon})$, $\nabla p \in W^{1,2}(B_{\varepsilon})$.

The second main result concerns the construction of the asymptotic approximations of the solution of problem (1.2). Let us describe the algorithm.

First, let us recall the definition of a quasi-Poiseuille flow for equations (1.2). Let $\sigma$ be a bounded domain with Lipschitz boundary in $\mathbb{R}^{n-1}$. Consider in the infinite cylinder $\Pi = \mathbb{R} \times \sigma$ the Dirichlet boundary value problem:

$$ \left\{ \begin{array}{l} -\mathrm{div}((\nu_0 + \lambda\nu(\dot{\gamma}(\mathbf{u})))D(\mathbf{u})) + \nabla p = 0, \quad x \in \Pi, \\ \mathrm{div}\mathbf{u} = 0, \quad x \in \Pi, \\ \mathbf{u}|_{\partial\Pi} = 0, \end{array} \right. \quad (1.4) $$

where $\dot{\gamma}(\mathbf{v}) = (d_{12}, d_{13}, d_{23}, 0, 0, 0)$ if $n = 3$ and $\dot{\gamma}(\mathbf{v}) = (d_{12}, 0, 0)$ if $n = 2$ (below we will see that for the quasi-Poiseuille flow $d_{ii} = 0$).

Define a quasi-Poiseuille flow as a solution to the following problem: find the couple $(\mathbf{V}_{P_{\alpha}}, \mathcal{P}_{P_{\alpha}})$ such that $\mathbf{V}_{P_{\alpha}}(x) = (v_{P_{\alpha}}(x'), 0, \dots, 0)^T$, and $\mathcal{P}_{P_{\alpha}}(x) = -\alpha x_1 + \beta$, $\alpha, \beta \in \mathbb{R}$, $x' = (x_2, \dots, x_n)$, where $v_{P_{\alpha}}$ is the solution of the following problem

$$ \left\{ \begin{array}{l} -\frac{1}{2} \mathrm{div}_{x'}((\nu_0 + \lambda \nu(\dot{\gamma}_P(v_{P_\alpha})))\nabla_{x'} v_{P_\alpha}) = \alpha, \quad x' \in \sigma, \\ v_{P_\alpha}|_{\partial \sigma} = 0. \end{array} \right. \quad (1.5) $$

Here $\dot{\gamma}_P(v_{P_\alpha}) = (\frac{1}{2}\nabla_{x'}v_{P_\alpha}, 0, 0)$ if $n = 2$, $\dot{\gamma}_P(v_{P_\alpha}) = (\frac{1}{2}\nabla_{x'}v_{P_\alpha}, 0, 0, 0)$ if $n = 3$, and $\alpha$ is the given pressure slope.

Define $F_{\sigma}(\alpha) = \int_{\sigma} v_{P_{\alpha}}(x')dx'$ the flux corresponding to the pressure slope $-\alpha$. Note that in the case of the steady Newtonian flow (the steady form of Navier-Stokes or Stokes equations) $F_{\sigma}(\alpha)$ is proportional to $\alpha$. This

<footnote n="1">¹In [27] the definition of $\dot{\gamma}(\mathbf{v})$ was introduced with forgotten $n$ last components. However it should be read as above.</footnote>
case corresponds to the value of $\lambda = 0$ and so, $F_\sigma(\alpha) = \kappa\alpha$, where $\kappa = \int_\sigma \tilde{v}_P(x')dx'$ and $\tilde{v}_P$ is a solution of the Poisson equation

$$ \left\{ \begin{array}{l} -\frac{\nu_0}{2} \Delta_{x'} \tilde{v}_P = 1, \quad x' \in \sigma, \\ \tilde{v}_P = 0, \quad x' \in \partial\sigma. \end{array} \right. \qquad (1.6) $$

Consider the quasi-Poiseuille flow in a thin tube with the cross-section $\sigma_\epsilon$ which is a contraction of $\sigma$ $1/\epsilon$-times. As the inflows/outflows have the velocity of order $\epsilon$ in $L^\infty$-norm, see (1.2)$_3$, we consider the Poiseuille velocity of the same order and it corresponds to the pressure slope of order $1/\epsilon$. So, for any $\alpha \in \mathbb{R}$ we consider the scaled problem (1.5) in $z' = \epsilon x'$ variables:

$$ \left\{ \begin{array}{l} -\frac{1}{2} \mathrm{div} z' \left( (\nu_0 + \lambda \nu (\dot{\gamma}_P(v_{P,\frac{\alpha}{\epsilon}}^\epsilon))) \nabla z' v_{P,\frac{\alpha}{\epsilon}}^\epsilon \right) = \frac{\alpha}{\epsilon}, \quad z' \in \sigma_\epsilon, \\ v_{P,\frac{\alpha}{\epsilon}}^\epsilon = 0, \quad z' \in \partial \sigma_\epsilon. \end{array} \right. \qquad (1.7) $$

Evidently, the solution of this problem is related to the solution of problem (1.5) as follows:

$$ v_{P,\frac{\alpha}{\epsilon}}^{\epsilon}(z') = \epsilon v_{P\alpha}(\frac{z'}{\epsilon}). \qquad (1.8) $$

Denote the corresponding flux

$$ F_{\sigma_{\epsilon}}(\frac{\alpha}{\epsilon}) = \int_{\sigma_{\epsilon}} v_{P,\frac{\alpha}{\epsilon}}^{\epsilon}(z') dz' = \int_{\sigma_{\epsilon}} \epsilon v_{P\alpha}(\frac{z'}{\epsilon}) dz' $$

After the change of variables $x' = \frac{z'}{\epsilon}$ we see that

$$ F_{\sigma_{\epsilon}}(\frac{\alpha}{\epsilon}) = \epsilon^n \int_{\sigma} v_{P\alpha}(x')dx' = \epsilon^n F_{\sigma}(\alpha), \qquad (1.9) $$

i.e.

$$ F_{\sigma_{\epsilon}}(\beta) = \epsilon^n F_{\sigma}(\epsilon\beta). \qquad (1.10) $$

These two last formulae give the scaling rule for the operator (function in stationary case) relating the pressure slope and the flux.<fnref n="2" />

In the same way one can introduce $\kappa_\epsilon = \epsilon^{n+1}\kappa$ the Newtonian flux in a contracted $1/\epsilon$-times cylinder corresponding to the pressure slope 1. Putting $G_{\sigma_\epsilon}(\beta) = F_{\sigma_\epsilon}(\beta) - \kappa_\epsilon\beta$, we get the relation

$$ G_{\sigma_{\epsilon}}(\beta) = \epsilon^n G_{\sigma}(\epsilon\beta). \qquad (1.11) $$

Consider now the problem on the graph corresponding to the data of problem (1.2): the cross-sections $\sigma_\epsilon^j$ and the given fluxes

<footnote n="2">²Note that the pressure slope in these problems is $-\alpha$ (or $-\alpha/\epsilon$) but in the name of the operator we skip the sign “minus”.</footnote>
$$F_l^\varepsilon = \int_{\gamma_l^\varepsilon} \varepsilon \mathbf{g}^l \left( \frac{x-O_l}{\varepsilon} \right) \cdot \mathbf{n}(x) dS = \varepsilon^n F^l,$$

$$ \left\{ \begin{array}{l} -\frac{\partial}{\partial x_1^{(e_j)}} \left( F_{\sigma_\varepsilon}^j \left( \frac{\partial p^\varepsilon}{\partial x_1^{(e_j)}}(x_1^{(e_j)}) \right) \right) = 0, \quad x_1^{(e_j)} \in (0, |e_j|), \\ -\sum_{e_j: O_l \in e_j} F_{\sigma_\varepsilon}^j \left( \frac{\partial p^\varepsilon}{\partial x_1^{(e_j)}}(0) \right) = 0, \quad l = 1, \dots, N_1, \\ -F_{\sigma_\varepsilon}^j \left( \frac{\partial p^\varepsilon}{\partial x_1^{(e_j)}}(0) \right) = -F_l^\varepsilon, \quad l = N_1 + 1, \dots, N, \quad O_l \in e_j, \\ p^\varepsilon(x_1^{(e_j)} = 0) - p^\varepsilon(x_1^{(e_s)} = 0) = 0, \quad e_j: O_l \in e_j, \quad l = 1, \dots, N_1, \\ p^\varepsilon(O_N) = 0. \end{array} \right. \qquad (1.12) $$

where $e_s$ is one of the edges with an end point $O_l$ called the selected edge of the node $O_l$. It can be rescaled using the scalings of the pressure slope - flux relation:

$$ \left\{ \begin{array}{l} -\frac{\partial}{\partial x_1^{(e_j)}} \left( F_{\sigma_j}(\varepsilon \frac{\partial p^\varepsilon}{\partial x_1^{(e_j)}}(x_1^{(e_j)})) \right) = 0, \quad x_1^{(e_j)} \in (0, |e_j|), \\ -\sum_{e_j: O_l \in e_j} F_{\sigma_j}(\varepsilon \frac{\partial p^\varepsilon}{\partial x_1^{(e_j)}}(0)) = 0, \quad l = 1, \dots, N_1, \\ -F_{\sigma_j}(\varepsilon \frac{\partial p^\varepsilon}{\partial x_1^{(e_j)}}(0)) = -F^l, \quad l = N_1 + 1, \dots, N, \quad O_l \in e_j, \\ \varepsilon p^\varepsilon(x_1^{(e_j)} = 0) - \varepsilon p^\varepsilon(x_1^{(e_s)} = 0) = 0, \quad e_j: O_l \in e_j, \quad l = 1, \dots, N_1, \\ \varepsilon p^\varepsilon(O_N) = 0. \end{array} \right. \qquad (1.13) $$

Notice that the operator $F_{\sigma_j}(\beta)$ is a nonlinear operator relating the pressure slope and the flux. Now making the change $p = \varepsilon p^\varepsilon$, we get a problem which does not depend on $\varepsilon$. If $p$ is an affine function, then

$$ p^{\varepsilon}(x_1^{(e_j)}) = -s_j x_1^{(e_j)} / \varepsilon + a_j / \varepsilon. \qquad (1.14) $$

Let us describe the leading term of the asymptotic expansion of the solution of problem (1.2). The pressure slope $\alpha_j^\varepsilon = -\frac{\partial p^\varepsilon}{\partial x_1^{(e_j)}}$ in every edge $e_j$ generates the quasi-Poiseuille velocity in the associated cylinder which is $v_{P,\alpha_j^\varepsilon}^\varepsilon(x^{(e_j)\prime})$. The functions $v_{P,\alpha_j^\varepsilon}^\varepsilon(x^{(e_j)\prime})$ (velocity) and $p^\varepsilon(x_1^{(e_j)}) - p^\varepsilon(x_1^{(e_s)} = 0)$ (difference of the pressure value at $x_1^{(e_j)}$ and of the pressure at the node $O_l$ for the end of the selected edge) are multiplied by a cut-off function $\zeta(\frac{x_1^{(e_j)}}{3r\varepsilon})$, where $\zeta$ is a $C^2$ smooth function equal to zero in the interval $[0, 1]$ and equal to one in the interval $[2, \infty)$. So, the regular part of the leading term is the couple $v_0^a, p_0^a$:

$$ \begin{array}{l} \mathbf{v}_0^a = v_{P,\alpha_j^\varepsilon}^\varepsilon(x^{(e_j)\prime}) \zeta(\frac{x_1^{(e_j)}}{3r\varepsilon}) e_j; \\ p_0^a = (p^\varepsilon(x_1^{(e_j)}) - p^\varepsilon(x_1^{(e_s)} = 0)) \zeta(\frac{x_1^{(e_j)}}{3r\varepsilon}) + p^\varepsilon(x_1^{(e_s)} = 0), \\ x_1^{(e_j)} \in (0, |e_j|/2). \end{array} \qquad (1.15) $$

where $e_j$ is the director vector of the edge $e_j$. This formula holds in the half of the edge $e_j$. For the other half we may use the analogous formula associated to the second end of the edge $e_j$.

This regular part of the leading term is completed with the boundary layer corrector also considered in the part $B_\varepsilon$ corresponding to $x_1^{(e_j)} \in (0, |e_j|/2)$. We say that the function $q \in L_{loc}^2(\Omega)$ exponentially stabilizes to
constants $q_1, q_2, \dots, q_J$ at infinity if

$$ \int_{\Omega^j} \exp(2\beta x_1^{(j)}) |q(x) - q_j|^2 dx < \infty, \quad j = 1, 2, \dots, J, $$

for some $\beta > 0$. The space of such functions we denote $\tilde{L}_{loc}^2(\Omega)$. Let $\mathbf{N}_0 \in \tilde{W}^{1,2}(\Omega_l)$ and $P_0 \in \tilde{L}_{loc}^2(\Omega_l)$ be solution of the following problem in dilated variables $\xi = \frac{x-O_l}{\epsilon}$ in the unbounded domain $\Omega_l$:

$$ \left\{ \begin{array}{l} -\operatorname{div}_{\xi}((\nu_0 + \lambda\nu(\dot{\gamma}_{\xi}(\mathbf{N}_0(\xi) + \mathbf{V}_{\xi}^0(\xi))))D_{\xi}(\mathbf{N}_0(\xi) + \mathbf{V}_{\xi}^0(\xi))) \\ \qquad + \nabla_{\xi}(P_0(\xi) + P_{\xi}^0(\xi)) = 0, \quad \xi \in \Omega_l, \\ \\ \operatorname{div}_{\xi}(\mathbf{N}_0(\xi) + \mathbf{V}_{\xi}^0(\xi)) = 0, \quad \xi \in \Omega_l, \\ \\ \mathbf{N}_0(\xi) = 0, \quad \xi \in \partial\Omega_l, \end{array} \right. \quad (1.16) $$

where

$$ \mathbf{V}_{\xi}^{0}(\xi) = \sum_{j:O_l \in e_j} \zeta(\frac{\xi_1^{(e_j)}}{3r})v_{P,s_j}(\xi^{(e_j)})e_j, \quad P_{\xi}^{0}(\xi) = \sum_{j:O_l \in e_j} \zeta(\frac{\xi_1^{(e_j)}}{3r})s_j\xi_1^{(e_j)}, $$

$\operatorname{div}_{\xi}, \dot{\gamma}_{\xi}, D_{\xi}, \nabla_{\xi}$ are operators written for $\xi$ variable. Notice that this problem is independent of $\epsilon$. Here and below for vertices $O_l$ condition $\mathbf{N}_k = 0$ on $\gamma^l$ is replaced by $\mathbf{N}_k = \mathbf{g}^l$.

So, the leading term with the boundary layer corrector has the form:

$$ \begin{aligned} \tilde{\mathbf{v}}_0^a &= v_{P,\alpha_j^\epsilon}^\epsilon(x^{(e_j)})\zeta(\frac{x_1^{(e_j)}}{3r\epsilon})e_j + \epsilon\mathbf{N}_0(\frac{x-O_l}{\epsilon}); \\ \tilde{p}_0^a &= (p^\epsilon(x_1^{(e_j)}) - p^\epsilon(x_1^{(e_s)} = 0))\zeta(\frac{x_1^{(e_j)}}{3r\epsilon}) + p^\epsilon(x_1^{(e_s)} = 0) + P^0(\frac{x-O_l}{\epsilon}), \\ x_1^{(e_j)} &\in (0, |e_j|/2). \end{aligned} \quad (1.17) $$

For the construction of high order asymptotic approximations of the solution we use a non-standard approach: instead of construction of asymptotic series as in [21] or [24] we construct a chain of successive iterations. Namely, in [26] it is proved that the pressure $P_0$ in (1.16) tends to some constants in the outlets. Taking into consideration that the pressure is defined up to an additive constant, we can assume that $P_0$ tends to zero at the outlet $\Pi_{\infty}^{(e_s)}$ corresponding to the selected edge $e_s$ of the bundle $B_l$. Denote by $\tilde{c}_{lj}^0$ the constants which are the limits of $P_0$ in the outlets $\Pi_{\infty}^{(e_j)}$, $j \neq s$. Then we define the first iteration of the problem on the graph (1.12) where the condition

$$ p^{\epsilon}(x_1^{(e_j)} = 0) - p^{\epsilon}(x_1^{(e_s)} = 0) = 0, \quad e_j: O_l \in e_j, \quad l = 1, \dots, N_1 $$

is replaced by

$$ p^{\epsilon}(x_1^{(e_j)} = 0) - p^{\epsilon}(x_1^{(e_s)} = 0) = \tilde{c}_{lj}^0, \quad e_j: O_l \in e_j, \quad l = 1, \dots, N_1. \quad (1.18) $$

The first iteration of the boundary layer corrector is the solution of problem (1.16) where in the expressions $\mathbf{V}_{\xi}^0$ and $P_{\xi}^0$ the constant $s_j$ is the slope of the solution of the problem on the graph $p^{\epsilon}$ of the first iteration. In turn, the pressure of this boundary layer problem tends to constants $\tilde{c}_{lj}^1$ and these constants appear as the right-hand side in of(1.18). This iterative procedure is continued. For the J-th iteration we get the estimate
of the error of order $O(\varepsilon^{J-2}|\ln \varepsilon|^{2J+2})$ in $W^{1,2}(B_\varepsilon)$-norm for the velocity and of order $O(\varepsilon^{J-3}|\ln \varepsilon|^{2J+2})$ in $L^2(B_\varepsilon)$-norm for the pressure.

The structure of the paper is as follows. In Section 2 the auxiliary domains are introduced: a covering of the tube structure and domains with cylindrical outlets to infinity. In Section 3 several embedding inequalities are proved for thin tube structures as well as a priori estimates for solutions of the divergence equation and Stokes problem posed in thin structure. Section 3 contains a generalized formulation of the Banach fixed point theorem used further for improving the regularity of solutions. Section 4 is devoted to the proof of the theorem on existence and uniqueness of a solution of the main problem for non-Newtonian flow and of a priori estimates for this solution. Sections 5-7 recall the results obtained in [25] and [26] on the steady non-Newtonian Poiseuille flow, on the problem on the graph, on the boundary layer problem in domains with cylindrical outlets. Finally, Section 8 is devoted to the construction of high order asymptotic approximations of the solution and the proof of the error estimates.

## 2. DEFINITIONS OF AUXILIARY DOMAINS

### 2.1. Covering of the domain $B_\varepsilon$

Let us construct a covering of the domain $B_\varepsilon$. Take domains $A_{\varepsilon,k}^{(e_j)} = \{x \in \Pi_\varepsilon^{(e_j)} : x_n^{(e_j)} \in \varepsilon(k-2, k+2)\}$, $j = 1, \dots, N, k = 2, \dots, L_\varepsilon^j, L_\varepsilon^j \sim |e|\varepsilon^{-1}$, and define $A_{\varepsilon,k}^{(j)} = \omega_\varepsilon^j \cup \{x \in \Pi_\varepsilon^{(e_j)} : x_n^{(e_j)} \in (0, 2\varepsilon)\}$, $j = N_1 + 1, \dots, N$ (i.e., when $O^j$ are vertices), and $A_{\varepsilon,k}^{(j)} = \omega_\varepsilon^j \cup \bigcup_{k_j} \{x \in \Pi_\varepsilon^{(e_{k_j})} : x_n^{(e_{k_j})} \in (0, 2\varepsilon)\}$, $j = 1, \dots, 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$. Obviously,

$$B_\varepsilon = \left( \bigcup_{j=1}^{N} \bigcup_{k=2}^{L_\varepsilon^j} A_{\varepsilon,k}^{(e_j)} \right) \cup \left( \bigcup_{j=1}^{N} A_{\varepsilon,k}^{(j)} \right).$$

We denote this covering by $\mathfrak{A}_\varepsilon$.

In parallel with the covering $\mathfrak{A}_\varepsilon$ we take the covering $\tilde{\mathfrak{A}}_\varepsilon$ containing larger domains

$$B_\varepsilon = \left( \bigcup_{j=1}^{N} \bigcup_{k=2}^{\tilde{L}_\varepsilon^j} \tilde{A}_{\varepsilon,k}^{(e_j)} \right) \cup \left( \bigcup_{j=1}^{N} \tilde{A}_{\varepsilon,k}^{(j)} \right),$$

where $\tilde{A}_{\varepsilon,k}^{(e_j)} = \{x \in \Pi_\varepsilon^{(e_j)} : x_n^{(e_j)} \in \varepsilon(k-3, k+3)\}$, $j = 1, \dots, N, k = 3, \dots, \tilde{L}_\varepsilon^j, \tilde{L}_\varepsilon^j \sim |e|\varepsilon^{-1}$. Then we define $A_{\varepsilon,k}^{(j)} = \omega_\varepsilon^j \cup \{x \in \Pi_\varepsilon^{(e_j)} : x_n^{(e_j)} \in (0, 3\varepsilon)\}$, $j = N_1 + 1, \dots, N$, and $A_{\varepsilon,k}^{(j)} = \omega_\varepsilon^j \cup \bigcup_{k_j} \{x \in \Pi_\varepsilon^{(e_{k_j})} : x_n^{(e_{k_j})} \in (0, 3\varepsilon)\}$, $j = 1, \dots, N_1$.

Obviously,

$$A_{\varepsilon,k}^{(e_j)} \subset \tilde{A}_{\varepsilon,k}^{(e_j)}, \quad A_{\varepsilon,k}^{(j)} \subset \tilde{A}_{\varepsilon,k}^{(j)}. \qquad (2.1)$$

### 2.2. Domains with cylindrical outlets to infinity

Consider the domain $\Omega \subset \mathbb{R}^n$ with $J$ cylindrical outlets to infinity which will be used in the construction of the boundary layer correctors for the asymptotic expansion of a solution of problem (1.2).

Let $\Omega = \Omega^0 \cup (\bigcup_{j=1}^J \Omega^j)$, where $\Omega^0$ is a bounded domain, $\Omega^0 \cap \Omega^j = \emptyset$ for $j \in \{1, \dots, J\}$, $\Omega^j \cap \Omega^l = \emptyset$ for $j \neq l$, $j, l \in \{1, \dots, J\}$, and the outlets to infinity $\Omega^j$ in some coordinate systems $x^{(j)} = (x_1^{(j)}, x^{(j)\prime})$, having the
origins within the boundary of domain $\Omega^0$, are given by the relations

$$ \Omega^j = \{x^{(j)} \in \mathbb{R}^n, x^{(j)\prime} \in \sigma_j, x_1^{(j)} \ge 0\}, $$

where $\sigma_j$ are some bounded domains in $\mathbb{R}^{n-1}$, cross-sections of the cylinders. Assume that for any $k \in \{1, \dots, J\}$ there exists a $\delta_j > 0$ such that the cylinder $\{x^{(j)} \in \mathbb{R}^n, x^{(j)\prime} \in \sigma_j, -\delta_j < x_1^{(j)} < 0\} \subset \Omega^0$. Denote $d_\sigma$ the maximal diameter of the cross-sections $\sigma_j$. We assume that the boundary $\partial\Omega$ is $C^4$-regular and that $\partial\Omega \cap \partial\Omega^0 \neq \emptyset$ has a positive measure. In particular, $\Omega$ can be just a semi-infinite cylinder: $\Omega = \{x \in \mathbb{R}^n, x' \in \sigma \subset \mathbb{R}^{n-1}, x_1 > 0\}$. Evidently there exists a positive real number $R > d_\sigma$ such that the ball $B_R = \{x \in \mathbb{R}^n, |x| < R\}$ contains $\Omega^0$. Note that $\Omega^j = \Pi_\infty^{(e)}$ for one of the edges of the bundle.

We also introduce the following notations:

$$ \Omega_{jk} = \{x \in \Omega_j : x_1^{(j)} < k\}, \quad \Omega^{(k)} = \Omega_0 \cup \left( \bigcup_{j=1}^{J} \Omega_{jk} \right), \qquad (2.2) $$

where $k \ge 0$ is an integer.

Let $\Omega \subset \mathbb{R}^n, n = 2, 3$, be domain with $J$ outlets to infinity. We define in $\Omega$ weighted function spaces. Denote $\beta = (\beta, \dots, \beta)$ and define by $\phi_\beta(x)$ a smooth function

$$ \phi_{\beta}(x) = \begin{cases} 0, & x \in \Omega_0, \\ \beta x_1^{(j)}, & x \in \Omega_j, x_1^{(j)} > 2, \quad j = 1, \dots, J. \end{cases} \qquad (2.3) $$

We also set $E_\beta(x) = \exp 2\phi_\beta(x)$.

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_1^{(j)} \to \infty$.

### 3. AUXILIARY RESULTS

#### 3.1. Embedding inequalities in tube structure $B_\epsilon$

**Lemma 3.1. (Poincaré inequality) There the inequality holds**

$$ \|u\|_{L^2(B_\epsilon)} \le c\epsilon \|\nabla u\|_{L^2(B_\epsilon)} \quad \forall u \in \dot{W}^{1,2}(B_\epsilon), \qquad (3.1) $$

where the constant $c$ is independent of $\epsilon$.

The proof of this lemma is obvious.

Let $\mathcal{G}$ be a bounded, Lipschitz domain in $\mathbb{R}^n$. Let us introduce in the Sobolev space $W^{l,2}(\mathcal{G})$ the equivalent (for the fixed $\epsilon$) norm

$$ \|u\|_{l,\alpha,\mathcal{G}}^2 = \sum_{k=0}^{l} \epsilon^{-2(l-k)+2\alpha} \|\nabla^k u\|_{L^2(\mathcal{G})}^2. $$
Obviously,

$$|||u|||_{l,\alpha,\mathcal{G}}^2 = \varepsilon^{2\alpha}|||u|||_{l,0,\mathcal{G}}^2.$$

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

$$\|u\|_{L^4(A_\varepsilon)}^4 \leq c\varepsilon^{-n+4} |||u|||_{1,0,A_\varepsilon}^4, \qquad (3.2)$$

where $A_{\varepsilon}$ is arbitrary domain from the covering $\mathcal{A}_{\varepsilon}$, and

$$\|u\|_{L^4(B_\varepsilon)}^4 \leq c\varepsilon^{-n+4} |||u|||_{1,0,B_\varepsilon}^4. \qquad (3.3)$$

*Proof.* In any bounded Lipschitz domain $\mathcal{G}$ the inequality holds (see [13])

$$\|u\|_{L^4(\mathcal{G})}^4 \leq c(\mathcal{G}) \left( \|u\|_{L^2(\mathcal{G})}^2 + \|\nabla u\|_{L^2(\mathcal{G})}^2 \right)^2.$$

By scaling, it is easy to see that in any $A_{\varepsilon}$ the estimate holds

$$\|u\|_{L^4(A_\varepsilon)}^4 \leq c\varepsilon^{-n+4} \left( \varepsilon^{-2}\|u\|_{L^2(A_\varepsilon)}^2 + \|\nabla u\|_{L^2(A_\varepsilon)}^2 \right)^2 = c\varepsilon^{-n+4} |||u|||_{1,0,A_\varepsilon}^4$$

with the constant $c$ independent of $\varepsilon$. From this it follows that

$$\|u\|_{L^4(A_\varepsilon)}^4 \leq c\varepsilon^{-n+4} |||u|||_{1,0,B_\varepsilon}^2 |||u|||_{1,0,A_\varepsilon}^2$$

and summing the last inequalities over all domains $A_{\varepsilon}$ from the covering $\mathcal{A}_{\varepsilon}$ we get (3.2). $\square$

By same token using the continuous embedding $W^{2,2}(\Omega)$ into $L^{\infty}(\Omega)$ we get the following
**Lemma 3.3.** Let $B_{\varepsilon} \subset \mathbb{R}^n$, $u \in W^{2,2}(B_{\varepsilon})$. Then the following inequalities

$$\begin{aligned} \|u\|_{L^{\infty}(A_{\varepsilon})} &\leq c\varepsilon^{-(n-4)/2} |||u|||_{2,0,A_{\varepsilon}}, \\ \|u\|_{L^{\infty}(B_{\varepsilon})} &\leq c\varepsilon^{-(n-4)/2} |||u|||_{2,0,B_{\varepsilon}} \end{aligned} \qquad (3.4)$$

hold with the constant $c$ independent of $\varepsilon$.

*Proof.* In a bounded Lipschitz domain $\mathcal{G}$ the inequality holds (see [13])

$$\|u\|_{L^{\infty}(\mathcal{G})}^2 \leq c(\mathcal{G}) \left( \|u\|_{L^2(\mathcal{G})}^2 + \|\nabla u\|_{L^2(\mathcal{G})}^2 + \|\nabla^2 u\|_{L^2(\mathcal{G})}^2 \right).$$

Then, by scaling, we obtain the estimate

$$\|u\|_{L^{\infty}(A_{\varepsilon})}^2 \leq c\varepsilon^{-n+4} |||u|||_{2,0,A_{\varepsilon}}^2$$

with the constant $c$ independent of $\varepsilon$. From this it follows that

$$\|u\|_{L^{\infty}(B_{\varepsilon})} \leq \sup_{A_{\varepsilon} \in \mathcal{A}_{\varepsilon}} \|u\|_{L^{\infty}(A_{\varepsilon})} \leq c\varepsilon^{-(n-4)/2} |||u|||_{2,0,B_{\varepsilon}}.$$

From Lemmas 3.2 and 3.3 follows $\square$
**Lemma 3.4.** Let $u \in W^{3,2}(B_\varepsilon)$. Then

$$ \|\nabla u\|_{L^\infty(B_\varepsilon)} \le c\varepsilon^{-(n-4)/2} \|u\|_{3,0, B_\varepsilon}, \quad (3.5) $$

$$ \|\nabla u\|_{L^4(B_\varepsilon)} \le c\varepsilon^{-(n-8)/4} \|u\|_{3,0, B_\varepsilon}, \quad (3.6) $$

$$ \|\nabla^2 u\|_{L^4(B_\varepsilon)} \le c\varepsilon^{-(n-4)/4} \|u\|_{3,0, B_\varepsilon}. \quad (3.7) $$

The same inequalities hold for functions $u \in W^{3,2}(A_\varepsilon)$.

## 3.2. Divergence equation

Let $G \subset \mathbb{R}^n$ be a bounded domain with Lipschitz boundary. Consider in $G$ the following problem:

For given $h \in L^2(G)$ with $\int_G h(x)dx = 0$ find a vector field $w \in \dot{W}^{1,2}(G)$ satisfying the equation

$$ \mathrm{div} w = h \quad \text{in} \quad G, \quad (3.8) $$

and the estimate

$$ \|\nabla w\|_{L^2(G)} \le c\|h\|_{L^2(G)}. \quad (3.9) $$

**Lemma 3.5.** Problem (3.8), (3.9) admits a solution. The constant $c$ in (3.9) depends only on the domain $G$.

Lemma 3.5 is proved in [12].

Consider now problem (3.8), (3.9) in the tube structure $B_\varepsilon$. The following result is obtained in [22].

**Lemma 3.6.** There exists a solution $w \in \dot{W}^{1,2}(B_\varepsilon)$ of problem (3.8), (3.9) in $B_\varepsilon$. There the estimate holds

$$ \|\nabla w\|_{L^2(B_\varepsilon)} \le c\varepsilon^{-1}\|h\|_{L^2(B_\varepsilon)}. \quad (3.10) $$

with the constant $c$ independent of $\varepsilon$.

Let us assume that $\partial B_\varepsilon \in C^4$ and introduce the notation $\hat{\Upsilon}_\varepsilon^{(e_j)} = \Pi_\varepsilon^{(e_j)} \cap \{x^{(e_j)} : x_1^{(e_j)} \in (\frac{1}{6}|e_j|, \frac{1}{3}|e_j|)\}$, $\hat{\Upsilon}_\varepsilon^{(e_j)} = \Pi_\varepsilon^{(e_j)} \cap \{x^{(e_j)} : x_1^{(e_j)} \in (\frac{1}{8}|e_j|, \frac{1}{2}|e_j|)\}$.

**Lemma 3.7.** Let $h \in W^{2,2}(\hat{\Upsilon}_\varepsilon^{(e_j)})$, $\mathrm{supp} h \subset \hat{\Upsilon}_\varepsilon^{(e_j)}$ and $\int_{\hat{\Upsilon}_\varepsilon^{(e_j)}} h(x)dx = 0$. Then there exists a solution $w \in \dot{W}^{1,2}(\hat{\Upsilon}_\varepsilon^{(e_j)}) \cap W^{3,2}(\hat{\Upsilon}_\varepsilon^{(e_j)})$ of the divergence equation (3.8). There the estimate holds

$$ \|w\|_{W^{3,2}(\hat{\Upsilon}_j)} \le c\varepsilon^{-3}\|h\|_{W^{2,2}(\hat{\Upsilon}_j)} \quad (3.11) $$

with the constant $c$ independent of $\varepsilon$.

Moreover, $\mathrm{supp} w \subset \{x : x_1 \in (\frac{1}{7}|e_j|, \frac{2}{5}|e_j|), x' \in \sigma^j\}$.

*Proof.* First consider the divergence equation in the domain $\hat{\Xi}_j = \{y : y_1 \in (\frac{1}{8}|e_j|, \frac{1}{2}|e_j|), y' \in \sigma^j\}$ assuming that $h \in W^{2,2}(\hat{\Xi}_j)$, $\mathrm{supp} h \subset \Xi_j = \{y : y_1 \in (\frac{1}{6}|e_j|, \frac{1}{3}|e_j|), y' \in \sigma^j\}$ and $\int_{\Xi_j} h(x)dx = 0$.
The lateral boundary of $\tilde{\Xi}_j$ is $C^4$-regular, $\text{supp} h \subset \Xi_j$, $\Xi_j \subset \tilde{\Xi}_j$ and $\Xi_j \neq \tilde{\Xi}_j$. So, by results in [8], there exits a vector field $\tilde{\mathbf{w}} \in \tilde{W}^{1,2}(\tilde{\Xi}_j) \cap W^{3,2}(\tilde{\Xi}_j)$ such that $\text{div}_y \tilde{\mathbf{w}}(y) = \tilde{h}(y)$, and

$$ \|\tilde{\mathbf{w}}\|_{W^{3,2}(\tilde{\Xi}_j)} \leq c \|\tilde{h}\|_{W^{2,2}(\tilde{\Xi}_j)}. \quad (3.12) $$

Moreover, $\tilde{\mathbf{w}}(y)$ can be constructed such that

$$ \text{supp} \tilde{\mathbf{w}} \subset \overline{\{y : y_1 \in (\frac{1}{7}|e_j|, \frac{2}{5}|e_j|), y' \in \sigma^j\}}. $$

Here $\tilde{h}(y) = h(x)|_{x=\mathcal{X}(y)}$, and $\mathcal{X}(y)$ is defined by $x_1 = y_1$, $x' = \varepsilon y'$.

Define the vector field $\mathbf{w}$ with components $w_1(x) = \tilde{w}_1(y)|_{\mathcal{X}^{-1}(x)}$, $w_i(x) = \varepsilon \tilde{w}_i(y)|_{\mathcal{X}^{-1}(x)}$, $i = 2, \dots, n$. It is straightforward to show that $\text{div}_x \mathbf{w}(x) = h(x)$, $\mathbf{w}|_{\tilde{\mathcal{Y}}_\varepsilon^{(e_j)}} = 0$. Passing in (3.12) to coordinates $x$ we obtain the inequality

$$ \int_{\tilde{\mathcal{Y}}_\varepsilon^{(e_j)}} \left( \sum_{l=0}^{3} \varepsilon^{2l} |\nabla_x^{\prime\prime} \tilde{\mathbf{w}}|^2 + \sum_{l=0}^{3} |\frac{\partial^l \tilde{\mathbf{w}}}{\partial x_1^l}|^2 \right) dx \leq c \int_{\tilde{\mathcal{Y}}_\varepsilon^{(e_j)}} \left( \sum_{l=0}^{2} \varepsilon^{2l} |\nabla_x^{\prime\prime} h|^2 + \sum_{l=0}^{2} |\frac{\partial^l h}{\partial x_1^l}|^2 \right) dx $$

which implies

$$ \sum_{l=0}^{3} \left( \varepsilon^{l-1} \|\nabla^l \mathbf{w}'\|_{L^2(\tilde{\mathcal{Y}}_\varepsilon^{(e_j)})} + \varepsilon^l \|\nabla^l w_1\|_{L^2(\tilde{\mathcal{Y}}_\varepsilon^{(e_j)})} + \varepsilon^{-1} \|\frac{\partial^l \mathbf{w}'}{\partial x_1^l}\|_{L^2(\tilde{\mathcal{Y}}_\varepsilon^{(e_j)})} + \|\frac{\partial^l w_1}{\partial x_1^l}\|_{L^2(\tilde{\mathcal{Y}}_\varepsilon^{(e_j)})} \right) \leq c \sum_{l=0}^{3} \left( \varepsilon^l \|\nabla^l h\|_{L^2(\tilde{\mathcal{Y}}_\varepsilon^{(e_j)})} + \|\frac{\partial^l h}{\partial x_1^l}\|_{L^2(\tilde{\mathcal{Y}}_\varepsilon^{(e_j)})} \right), $$

from which we obtain (3.12). $\square$

### 3.3. Stokes problem

Denote by $H(B_\varepsilon)$ the subspace of divergence free functions from $\tilde{W}^{1,2}(B_\varepsilon)$.

Consider in $B_\varepsilon$ the Dirichlet problem for the Stokes system

$$ \left\{ \begin{array}{ll} -\nu \Delta \mathbf{v} + \nabla p = \mathbf{f}, & x \in B_\varepsilon, \\ \text{div} \mathbf{v} = 0, & x \in B_\varepsilon, \\ \mathbf{v}|_{\partial B_\varepsilon} = 0. & \end{array} \right. \qquad (3.13) $$

The weak solution $\mathbf{v} \in H(B_\varepsilon)$ to (3.13) satisfies the integral identity

$$ \nu \int_{B_\varepsilon} \nabla \mathbf{v} \cdot \nabla \eta \, dx = \int_{B_\varepsilon} \mathbf{f} \cdot \eta \, dx \quad \forall \eta \in H(B_\varepsilon), $$

and hence the estimate

$$ \|\nabla \mathbf{v}\|_{L^2(B_\varepsilon)}^2 \leq c \varepsilon^2 \|\mathbf{f}\|_{L^2(B_\varepsilon)}^2. \qquad (3.14) $$

**Lemma 3.8.** Let $\partial B_\varepsilon \in C^{l+2}$, $\mathbf{f} \in W^{l,2}(B_\varepsilon)$. Then $\mathbf{v} \in W^{l+2,2}(B_\varepsilon)$, $\nabla p \in W^{l,2}(B_\varepsilon)$ and

$$ \|\mathbf{v}\|_{l+2,\alpha, B_\varepsilon}^2 + \|\nabla p\|_{l,\alpha, B_\varepsilon}^2 \leq c \|\mathbf{f}\|_{l,\alpha, B_\varepsilon}^2. \qquad (3.15) $$
with the constant c independent of $\varepsilon$.

*Proof.* Let $A_\varepsilon \subset \tilde{A}_\varepsilon$ be domains from the covering $\mathfrak{A}_\varepsilon$ and $\tilde{\mathfrak{A}}_\varepsilon$ of $B_\varepsilon$. Consider (3.13) in $\tilde{A}_\varepsilon$. Making the change of variables $y = \varepsilon^{-1}x$ we transform $A_\varepsilon$ and $\tilde{A}_\varepsilon$ into the fixed (independent of $\varepsilon$) domains $A_0$ and $\tilde{A}_0$. The Stokes problem in coordinates $y$ takes the form

$$ egin{cases} -\nu \Delta_y \mathbf{v} + \nabla_y (\varepsilon p) = \varepsilon^2 \mathbf{f} & \text{in } \tilde{A}_0 \\ \operatorname{div}_y \mathbf{v} = 0 & \text{in } \tilde{A}_0, \\ \mathbf{v}|_{\partial B_\varepsilon \cap \partial \tilde{A}_0} = 0. \end{cases} \tag{3.16} $$

ADN local estimates for elliptic problems (see [2]) yield the inequality

$$ \sum_{m=0}^{l+2} \|\nabla^m \mathbf{v}\|_{L^2(\tilde{A}_0)}^2 + \sum_{m=1}^{l+1} \|\nabla^m q\|_{L^2(\tilde{A}_0)}^2 \\ \le c \left( \varepsilon^4 \sum_{m=0}^{l} \|\nabla^m \mathbf{f}\|_{L^2(\tilde{A}_0)}^2 + \|\mathbf{v}\|_{L^2(\tilde{A}_0)}^2 + \|q - \bar{q}\|_{L^2(\tilde{A}_0)}^2 \right), \tag{3.17} $$

where $q(y) = \varepsilon p(y)$, $\bar{q} = \frac{1}{|\tilde{A}_0|} \int_{\tilde{A}_0} q(y) dy$. Since $\int_{\tilde{A}_0} (q(y) - \bar{q}) dy = 0$, there exists $\mathbf{w} \in \dot{W}^{1,2}(\tilde{A}_0)$ such that $\operatorname{div} \mathbf{w} = q(y) - \bar{q}$ in $\tilde{A}_0$ and

$$ \|\nabla \mathbf{w}\|_{L^2(\tilde{A}_0)} \le c \|q - \bar{q}\|_{L^2(\tilde{A}_0)}. $$

Multiplying (3.16) by $\mathbf{w}$ and integrating by parts we obtain

$$ \begin{aligned} \|q - \bar{q}\|_{L^2(\tilde{A}_0)}^2 &= \int_{\tilde{A}_0} q(y)(q(y) - \bar{q})dy = \int_{\tilde{A}_0} q(y)\operatorname{div}\mathbf{w}dy \\ &= \nu \int_{\tilde{A}_0} \nabla\mathbf{v} \cdot \nabla\mathbf{w}dy - \varepsilon^2 \int_{\tilde{A}_0} \mathbf{f} \cdot \mathbf{w}dy \\ &\le \|\nabla\mathbf{v}\|_{L^2(\tilde{A}_0)} \|\nabla\mathbf{w}\|_{L^2(\tilde{A}_0)} + \varepsilon^2 \|\mathbf{f}\|_{L^2(\tilde{A}_0)} \|\mathbf{w}\|_{L^2(\tilde{A}_0)} \\ &\le c \|\nabla\mathbf{v}\|_{L^2(\tilde{A}_0)} \|q - \bar{q}\|_{L^2(\tilde{A}_0)} + c\varepsilon^2 \|\mathbf{f}\|_{L^2(\tilde{A}_0)} \|q - \bar{q}\|_{L^2(\tilde{A}_0)}. \end{aligned} $$

Therefore,

$$ \|q - \bar{q}\|_{L^2(\tilde{A}_0)} \le c (\|\nabla \mathbf{v}\|_{L^2(\tilde{A}_0)} + \varepsilon^2 \|\mathbf{f}\|_{L^2(\tilde{A}_0)}). \tag{3.18} $$

From (3.17), using (3.18) and the Poincaré inequality, we derive

$$ \sum_{m=0}^{l+2} \|\nabla^m \mathbf{v}\|_{L^2(\tilde{A}_0)}^2 + \sum_{m=1}^{l+1} \|\nabla^m q\|_{L^2(\tilde{A}_0)}^2 \\ \le c \left( \varepsilon^4 \sum_{m=0}^{l} \|\nabla^m \mathbf{f}\|_{L^2(\tilde{A}_0)}^2 + \|\nabla_y \mathbf{v}\|_{L^2(\tilde{A}_0)}^2 \right). \tag{3.19} $$

Returning to coordinates $x$ we obtain

$$ \begin{aligned} & \sum_{m=0}^{l+2} \varepsilon^{2m} \|\nabla_x^m \mathbf{v}\|_{L^2(\tilde{A}_\varepsilon)}^2 + \varepsilon^2 \sum_{m=1}^{l+1} \varepsilon^{2m} \|\nabla_x^m p\|_{L^2(\tilde{A}_\varepsilon)}^2 \\ & \le c \left( \varepsilon^4 \sum_{m=0}^{l} \varepsilon^{2m} \|\mathbf{f}\|_{L^2(\tilde{A}_\varepsilon)}^2 + \varepsilon^2 \|\nabla \mathbf{v}\|_{L^2(\tilde{A}_\varepsilon)}^2 \right). \end{aligned} \tag{3.20} $$
Summing up (3.20) over all domains $A_{\varepsilon} \subset \tilde{A}_{\varepsilon}$ and estimating the last term in the right hand side by (3.14) yields

$$ \sum_{m=0}^{l+2} \varepsilon^{2m} \|\nabla_x^m \mathbf{v}\|_{L^2(B_\varepsilon)}^2 + \varepsilon^2 \sum_{m=1}^{l+1} \varepsilon^{2m} \|\nabla_x^m p\|_{L^2(B_\varepsilon)}^2 \\ \le c \left( \varepsilon^4 \sum_{m=0}^{l} \varepsilon^{2m} \|\mathbf{f}\|_{L^2(B_\varepsilon)}^2 + \varepsilon^2 \|\nabla_x \mathbf{v}\|_{L^2(B_\varepsilon)}^2 \right) \le c \varepsilon^4 \sum_{m=0}^{l} \varepsilon^{2m} \|\mathbf{f}\|_{L^2(B_\varepsilon)}^2. \quad (3.21) $$

Multiplying the last inequality by $\varepsilon^{-2(l+2)+2\alpha}$ we obtain (3.15). □

Consider now the Stokes problem in the domain $\Omega$ with $J$ outlets to infinity:

$$ \left\{ \begin{array}{ll} -\nu \Delta \mathbf{v} + \nabla p = \mathbf{f}, & x \in \Omega, \\ \operatorname{div} \mathbf{v} = 0, & x \in \Omega, \\ \mathbf{v}|_{\partial \Omega} = 0. \end{array} \right. \quad (3.22) $$

We have the weak solution $\mathbf{v} \in H(\Omega)$ to (3.22) which satisfies the integral identity

$$ \nu \int_{\Omega} \nabla \mathbf{v} \cdot \nabla \boldsymbol{\eta} \mathrm{d}x = \int_{\Omega} \mathbf{f} \cdot \boldsymbol{\eta} \mathrm{d}x \quad \forall \boldsymbol{\eta} \in H(\Omega), $$

and the estimate

$$ \|\nabla \mathbf{v}\|_{L^2(\Omega)}^2 \le c \|\mathbf{f}\|_{L^2(\Omega)}^2. \quad (3.23) $$

The following theorem was proved in [30] (see also [28], Thm. III.3.2).

**Theorem 3.9.** Let $\partial\Omega \in C^{l+2}$ and let $\mathbf{f} \in W_{2,\beta}^l(\Omega)$, $\beta \ge 0$. Suppose that the number $0 \le \beta \le \beta_*$, If $\beta_*$ is sufficiently small, then the weak solution $\mathbf{u}$ belongs to the space $W_{2,\beta}^{l+2}(\Omega)$ and there exists a pressure function $p$ with $\nabla p \in W_{2,\beta}^l(\Omega)$ such that the pair $(\mathbf{u}(x), p(x))$ satisfies equations (3.1) almost everywhere in $\Omega$. There the estimate holds

$$ \|\mathbf{u}\|_{W_{2,\beta}^{l+2}(\Omega)} + \|\nabla p\|_{W_{2,\beta}^l(\Omega)} \le c \|\mathbf{f}\|_{W_{2,\beta}^l(\Omega)}. \quad (3.24) $$

### 3.4. Weak Banach contraction principle

**Theorem 3.10.** Let $X$ and $Y$ be reflexive Banach spaces, $X \subset Y$,

$$ \|x\|_Y \le \|x\|_X \quad \forall x \in X. \quad (3.25) $$

Suppose that $M \subset X$ is closed, bounded set, $M \ne \emptyset$, and the mapping $T: M \to M$ satisfies the inequality

$$ \|Tx - Ty\|_Y \le k\|x - y\|_Y, \quad k < 1. \quad (3.26) $$

Then $T$ admits exactly one fixed point $x_* \in M$:

$$ Tx_* = x_* $$

This result is well known and widely used in the mathematical community. For the proof of it see, for example, [26].
4. EXISTENCE AND UNIQUENESS OF THE SOLUTION OF THE MAIN PROBLEM

Since $\mathrm{div} \mathbf{v} = 0$, the problem (1.2) can be written in the form

$$ egin{cases} -\frac{\nu_0}{2} \Delta \mathbf{v} + \nabla p = -\lambda \mathrm{div} \left( \nu(\dot{\gamma}(\mathbf{v})) D(\mathbf{v}) \right), \\ \mathrm{div} \mathbf{v} = 0, \quad x \in B_{\varepsilon} \\ \mathbf{v}|_{\partial B_{\varepsilon}} = \varepsilon \mathbf{g}, \end{cases} \quad (4.1) $$

Let $e = e_{O_j}$ be the edge with the end $O_j$ and let $x^{(e)}$ be the Cartesian coordinates corresponding to the origin $O_j$ and the edge $e$, i.e., $x^{(e)} = \mathcal{P}^{(e)}(x - O_j)$, $\mathcal{P}^{(e)}$ is the orthogonal matrix relating the global coordinates $x$ with the local ones $x^{(e)}$, $\sigma_{\varepsilon}^j = \{x : \frac{x^{(e)}}{\varepsilon} \in \sigma, x_1^{(e)} = 0\}$. Denote $\mathbf{g}^{(e)} = \mathcal{P}^{(e)} \mathbf{g}^j$.

Let

$$ egin{aligned} \tilde{F}^j &= \varepsilon \int \mathbf{g}(x) \cdot \mathbf{n}(x) dS = \varepsilon \int \mathbf{g}^j \left( \frac{x - O_j}{\varepsilon} \right) \cdot \mathbf{n}(x) dS \\ &= \varepsilon^n \int \hat{g}_n^j(y^{(e)}) dy^{(e)} \equiv \varepsilon^n F^j, \quad j = N_1 + 1, \dots, N, \end{aligned} \quad (4.2) $$

where $\mathbf{n}$ is the unit outward (with respect to $B_{\varepsilon}$) normal vector to $\gamma_{\varepsilon}^j$, $y^{(e)} = \frac{x^{(e)}}{\varepsilon}$, $\hat{g}^j(y^{(e)}) = \mathbf{g}^j((\mathcal{P}^{(e)})^*y^{(e)})$, $F^j$ does not depend on $\varepsilon$. Assume that for the flow rates $F^j$ the compatibility condition

$$ \sum_{j=N_1+1}^{N} F^j = 0 \quad (4.3) $$

is valid.

Let $\mathbf{g}$ be the divergence free extension of the boundary function $\mathbf{g}$ (which we denote by the same symbol $\mathbf{g}$, $\mathbf{g} \in W^{3,2}(B_{\varepsilon})$) satisfying the following asymptotic estimates

$$ \|\nabla^l \mathbf{g}\|_{L^2(A_\varepsilon)} \leq c G_0 \varepsilon^{\frac{n-2l}{2}}, \quad l = 0, 1, 2, 3, \quad (4.4) $$

for any domain $A_{\varepsilon}$ from the covering $\mathfrak{A}_{\varepsilon}$. Using Lemma 3.3 we see that

$$ \begin{gathered} \sup_{x \in B_{\varepsilon}} |\mathbf{g}(x)| \le cG_0, \quad \sup_{x \in B_{\varepsilon}} |\nabla \mathbf{g}(x)| \le cG_0 \varepsilon^{-1}, \\ \|\nabla^l \mathbf{g}\|_{L^2(B_{\varepsilon})} \le cG_0 \varepsilon^{\frac{n-(2l+1)}{2}}, \quad l = 0, 1, 2, 3, \end{gathered} \quad (4.5) $$

where the constant $c$ is independent of $\varepsilon$ and $\mathbf{g}$ and $G_0$ is independent of $\varepsilon$.

Representing $\mathbf{v}$ as the sum $\mathbf{v} = \mathbf{u} + \varepsilon \mathbf{g}$ we obtain the following problem

$$ egin{cases} -\frac{\nu_0}{2} \Delta \mathbf{u} + \nabla p = \lambda \mathrm{div} \left( \nu(\dot{\gamma}(\mathbf{u} + \varepsilon \mathbf{g})) D(\mathbf{u} + \varepsilon \mathbf{g}) \right) \\ -\lambda \mathrm{div} \left( \nu(\dot{\gamma}(\varepsilon \mathbf{g})) D(\varepsilon \mathbf{g}) \right) + \mathbf{f}, \\ \mathrm{div} \mathbf{u} = 0, \quad x \in B_{\varepsilon}, \\ \mathbf{u}|_{\partial B_{\varepsilon}} = 0, \end{cases} \quad (4.6) $$

where $\mathbf{f} = \frac{\varepsilon \nu_0}{2} \Delta \mathbf{g} + \lambda \mathrm{div} \left( \nu(\dot{\gamma}(\varepsilon \mathbf{g})) D(\varepsilon \mathbf{g}) \right)$. Below we consider problem (4.6) with arbitrary $\mathbf{f} \in W^{1,2}(B_{\varepsilon})$.
**Theorem 4.1.** Let **f** be a vector-valued function from $W^{1,2}(B_\varepsilon)$ and $\mathbf{g} \in W^{3,2}(B_\varepsilon)$ be a described above extension.

(i) There exists $\lambda_0$ such that for all $\lambda \in (0, \lambda_0)$ problem (4.6) admits a unique in some ball solution $(\mathbf{u}, p)$ with $\mathbf{u} \in W^{3,2}(B_\varepsilon)$, $\nabla p \in W^{1,2}(B_\varepsilon)$. There holds the estimate

$$ |||\mathbf{u}|||_{3,0,B_\varepsilon}^2 + |||\nabla p|||_{1,0,B_\varepsilon}^2 \le c |||\mathbf{f}|||_{1,0,B_\varepsilon}^2. \quad (4.7) $$

Here and below constants $c$ with or without subscripts are independent of $\varepsilon$.

(ii) The pressure $p$ is unique up to an additive constant. Being normalized by the condition of the zero mean value $\int_{B_\varepsilon} pdx = 0$, the pressure satisfies the estimate

$$ \|p\|_{L^2(B_\varepsilon)}^2 \le c\varepsilon^2 |||\mathbf{f}|||_{1,0,B_\varepsilon}^2. \quad (4.8) $$

*Proof.* (i) Let $\mathcal{M}$ be the operator $H(B_\varepsilon) \cap W^{3,2}(B_\varepsilon) \to H(B_\varepsilon) \cap W^{3,2}(B_\varepsilon)$, such that for any $\mathbf{U} \in H(B_\varepsilon) \cap W^{3,2}(B_\varepsilon)$, $(\mathcal{M}\mathbf{U}, P)$ is a solution of the problem

$$ \left\{ \begin{array}{ll} -\frac{\nu_0}{2} \Delta \mathcal{M} \mathbf{U} + \nabla P = \mathbf{h}(\mathbf{U} + \varepsilon \mathbf{g}) + \mathbf{f}, & x \in B_\varepsilon \\ \mathrm{div} \mathcal{M} \mathbf{U} = 0, & x \in B_\varepsilon \\ \mathcal{M} \mathbf{U}|_{\partial B_\varepsilon} = 0, & \end{array} \right. \quad (4.9) $$

where

$$ \begin{aligned} \mathbf{h}(\mathbf{U} + \varepsilon \mathbf{g}) &= \lambda \mathrm{div} \left( \nu(\dot{\gamma}(\mathbf{U} + \varepsilon \mathbf{g})) D(\mathbf{U} + \varepsilon \mathbf{g}) \right) - \lambda \mathrm{div} \left( \nu(\dot{\gamma}(\varepsilon \mathbf{g})) D(\varepsilon \mathbf{g}) \right) \\ &= \lambda \mathrm{div} \left[ \nu(\dot{\gamma}(\mathbf{U} + \varepsilon \mathbf{g})) D(\mathbf{U}) + \left( \nu(\dot{\gamma}(\mathbf{U} + \varepsilon \mathbf{g})) - \nu(\dot{\gamma}(\varepsilon \mathbf{g})) \right) D(\varepsilon \mathbf{g}) \right] \\ &= \lambda \nabla^T \nu(\dot{\gamma}(\mathbf{U} + \varepsilon \mathbf{g})) \cdot D(\mathbf{U}) + \lambda \nu(\dot{\gamma}(\mathbf{U} + \varepsilon \mathbf{g})) \mathrm{div} D(\mathbf{U}) \\ &\quad + \lambda \left( \nabla^T \nu(\dot{\gamma}(\mathbf{U} + \varepsilon \mathbf{g})) - \nabla^T \nu(\dot{\gamma}(\varepsilon \mathbf{g})) \right) \cdot D(\varepsilon \mathbf{g}) + \lambda \left( \nu(\dot{\gamma}(\mathbf{U} + \varepsilon \mathbf{g})) - \nu(\dot{\gamma}(\varepsilon \mathbf{g})) \right) \mathrm{div} D(\varepsilon \mathbf{g}). \end{aligned} $$

Here and below the gradient $\nabla$ is a column vector.

Note that

$$ \nabla \nu(\dot{\gamma}(\mathbf{U} + \varepsilon \mathbf{g})) = (\nabla_y \nu(y)|_{y=\dot{\gamma}(\mathbf{U}+\varepsilon \mathbf{g})})^T \nabla \dot{\gamma}(\mathbf{U} + \varepsilon \mathbf{g}) $$

and

$$ \nabla \nu(\dot{\gamma}(\varepsilon \mathbf{g})) = (\nabla_y \nu(y)|_{y=\dot{\gamma}(\varepsilon \mathbf{g})})^T \nabla \dot{\gamma}(\varepsilon \mathbf{g}), $$

where $\nabla \dot{\gamma}$ is the Jacobian matrix of $\dot{\gamma}$. Let us subtract and add in the term $\lambda(\nabla^T \nu(\dot{\gamma}(\mathbf{U} + \varepsilon \mathbf{g})) - \nabla^T \nu(\dot{\gamma}(\varepsilon \mathbf{g})))$ the expression

$$ \lambda(\nabla_y \nu(y)|_{y=\dot{\gamma}(\mathbf{U}+\varepsilon \mathbf{g})})^T \nabla^T \dot{\gamma}(\varepsilon \mathbf{g}) $$
and notice that

$$|\lambda(\nabla_y \nu(y)|_{y=\dot{\gamma}(\mathbf{U}+\varepsilon\mathbf{g})})^T \nabla^T \dot{\gamma}(\mathbf{U} + \varepsilon\mathbf{g})) - \lambda(\nabla_y \nu(y)|_{y=\dot{\gamma}(\mathbf{U}+\varepsilon\mathbf{g})})^T \nabla^T \dot{\gamma}(\varepsilon\mathbf{g})| \\ \leq c\lambda \sup_y |\nabla_y \nu(y)||\nabla^2 U|$$

and

$$|\lambda(\nabla_y \nu(y)|_{y=\dot{\gamma}(\mathbf{U}+\varepsilon\mathbf{g})})^T \nabla^T \dot{\gamma}(\varepsilon\mathbf{g})) - \lambda(\nabla_y \nu(y)|_{y=\dot{\gamma}(\varepsilon\mathbf{g})})^T \nabla^T \dot{\gamma}(\varepsilon\mathbf{g})| \\ \leq c\lambda \sup_y |\nabla_y^2 \nu(y)||\nabla \dot{\gamma}(\varepsilon\mathbf{g})||\nabla U|.$$

These estimates will be used below.

Let us estimate the right hand side of equation (4.9). By (1.1) and (4.5), we have

$$|\mathbf{h}| \leq c\lambda \sup_y |\nabla_y \nu(y)||\nabla \mathbf{U}||\nabla^2(\mathbf{U} + \varepsilon\mathbf{g})| + c\lambda \sup_y |\nu(y)||\nabla^2 \mathbf{U}|$$

$$+ c\lambda \sup_y |\nabla_y \nu(y)||\nabla \mathbf{U}||\varepsilon \nabla^2 \mathbf{g}| + c\lambda \sup_y |\nabla_y \nu(y)||\nabla^2 \mathbf{U}||\varepsilon \nabla \mathbf{g}|$$

$$+ c\lambda \sup_y |\nabla_y^2 \nu(y)||\nabla \mathbf{U}||\varepsilon \nabla^2 \mathbf{g}||\varepsilon \nabla \mathbf{g}|$$

$$\leq c\lambda A (|\nabla \mathbf{U}||\nabla^2 \mathbf{U}| + |\nabla \mathbf{U}||\varepsilon \nabla^2 \mathbf{g}| + |\nabla^2 \mathbf{U}| + |\nabla^2 \mathbf{U}||\varepsilon \nabla \mathbf{g}| + |\nabla \mathbf{U}||\varepsilon \nabla^2 \mathbf{g}||\varepsilon \nabla \mathbf{g}|).$$

Using (3.5) we obtain

$$\int_{B_{\varepsilon}} |\nabla^2 \mathbf{U}|^2 |\nabla \mathbf{U}|^2 dx \leq \sup_{x \in B_{\varepsilon}} |\nabla \mathbf{U}(x)|^2 \int_{B_{\varepsilon}} |\nabla^2 \mathbf{U}|^2 dx \leq c \varepsilon^{2\alpha_1+2} |||\mathbf{U}|||_{3,0, B_{\varepsilon}}^4,$$

where $\alpha_1 = 1$ for $n = 2$ and $\alpha_1 = 1/2$ for $n = 3$. Further, applying (4.5) we get

$$\varepsilon^2 \int_{B_{\varepsilon}} |\nabla \mathbf{g}|^2 |\nabla^2 \mathbf{U}|^2 dx \leq \sup_{x \in B_{\varepsilon}} |\nabla \mathbf{g}(x)|^2 \int_{B_{\varepsilon}} |\nabla^2 \mathbf{U}|^2 \leq c G_0^2 \varepsilon^2 |||\mathbf{U}|||_{3,0, B_{\varepsilon}}^2.$$

Let us estimate the integral containing the term $\varepsilon^2|\nabla^2 \mathbf{g}|^2|\nabla \mathbf{U}|^2$. First consider this integral over the domain $A_{\varepsilon}$, where $A_{\varepsilon}$ is an arbitrary domain from the covering $\mathfrak{A}_{\varepsilon}$. Inequalities (3.5), (4.4) yield

$$\varepsilon^2 \int_{A_{\varepsilon}} |\nabla^2 \mathbf{g}|^2 |\nabla \mathbf{U}|^2 dx \leq c \varepsilon^2 \sup_{x \in A_{\varepsilon}} |\nabla \mathbf{U}(x)|^2 \int_{A_{\varepsilon}} |\nabla^2 \mathbf{g}|^2 dx$$

$$\leq c \varepsilon^{2+2\alpha_1} |||\mathbf{U}|||_{3,0, A_{\varepsilon}}^2 G_0^2 \varepsilon^{n-4} \leq c G_0^2 \varepsilon^2 |||\mathbf{U}|||_{3,0, A_{\varepsilon}}^2;$$
$$ \varepsilon^4 \int_{A_\varepsilon} |\nabla \mathbf{g}|^2 |\nabla^2 \mathbf{g}|^2 |\nabla \mathbf{U}|^2 dx \le c \varepsilon^4 \sup_{x \in A_\varepsilon} |\nabla \mathbf{g}(x)|^2 \int_{A_\varepsilon} |\nabla^2 \mathbf{g}|^2 |\nabla \mathbf{U}|^2 dx $$

$$ \le c \varepsilon^4 \sup_{x \in B_\varepsilon} |\nabla \mathbf{g}(x)|^2 \int_{A_\varepsilon} |\nabla^2 \mathbf{g}|^2 |\nabla \mathbf{U}|^2 dx \le c G_0^4 \varepsilon^2 |||\mathbf{U}|||_{3,0,A_\varepsilon}^2 $$

Summing these inequalities over all $A_\varepsilon \in \mathfrak{A}_\varepsilon$, we derive

$$ \varepsilon^2 \int_{B_\varepsilon} |\nabla^2 \mathbf{g}|^2 |\nabla \mathbf{U}|^2 dx + \varepsilon^4 \int_{B_\varepsilon} |\nabla^2 \mathbf{g}|^2 |\nabla^2 \mathbf{g}|^2 |\nabla \mathbf{U}|^2 dx $$

$$ \le c G_0^2 (1 + G_0^2) \varepsilon^2 |||\mathbf{U}|||_{3,0,B_\varepsilon}^2 $$

Collecting the above estimates we get

$$ \frac{1}{\varepsilon^2} \|\mathbf{h}\|_{L^2(B_\varepsilon)}^2 \le c \lambda^2 A^2 \left( (1 + G_0^2(1 + G_0^2)) |||\mathbf{U}|||_{3,0,B_\varepsilon}^2 + \varepsilon^{2\alpha_1} |||\mathbf{U}|||_{3,0,B_\varepsilon}^4 \right). \quad (4.10) $$

Analogously, we have

$$ \begin{aligned} |\nabla \mathbf{h}| &\le c \lambda A \left( (|\nabla^3 \mathbf{U}| + |\nabla^2 \mathbf{U}|^2 + \varepsilon |\nabla^3 \mathbf{U}||\nabla \mathbf{g}|) (1 + |\nabla \mathbf{U}|) \right. \\ &\qquad \left. + (\varepsilon |\nabla^2 \mathbf{U}||\nabla^2 \mathbf{g}| + \varepsilon |\nabla^2 \mathbf{U}|^2 |\nabla \mathbf{g}|) (1 + |\nabla \mathbf{U}|) + \varepsilon^2 |\nabla^2 \mathbf{U}||\nabla^2 \mathbf{g}||\nabla \mathbf{g}| \right. \\ &\qquad \left. + |\nabla \mathbf{U}| (\varepsilon |\nabla^3 \mathbf{g}| + \varepsilon^2 |\nabla^2 \mathbf{g}|^2) (1 + \varepsilon |\nabla \mathbf{g}|) \right). \end{aligned} $$

The $L^2$ norm of this expression is evaluated according to the following scheme: in each product of gradients the first order terms $|\nabla \mathbf{U}|$ and $\varepsilon |\nabla \mathbf{g}|$ are evaluated by $\sup_{x \in B_\varepsilon} |\nabla \mathbf{U}(x)|$ and $\sup_{x \in B_\varepsilon} \varepsilon |\nabla \mathbf{g}|$, the second order terms $|\nabla^2 \mathbf{U}|$ and $\varepsilon |\nabla \mathbf{g}|$ are evaluated in the $L^4$ norm, finally the third order terms $|\nabla^3 \mathbf{U}|$ and $\varepsilon |\nabla \mathbf{g}|$ are evaluated in the $L^2$ norm. Then we apply the embedding inequalities of Lemma 3.4. So, for the gradient of $\mathbf{h}$ using (1.1), (4.4), (4.5), (3.5)-(3.7) we obtain the estimate<fnref n="3" />

$$ \|\nabla \mathbf{h}\|_{L^2(B_\varepsilon)}^2 \le c \lambda^2 A^2 (1 + G_0^2 + G_0^4) (|||\mathbf{U}|||_{3,0,B_\varepsilon}^6 + |||\mathbf{U}|||_{3,0,B_\varepsilon}^4 + |||\mathbf{U}|||_{3,0,B_\varepsilon}^2). \quad (4.11) $$

Let us define in $W^{3,2}(B_\varepsilon)$ a closed bounded set $B_{R_0} = \{\mathbf{u} \in W^{3,2}(B_\varepsilon) : |||\mathbf{u}|||_{3,0,B_\varepsilon} \le R_0\}$. Assume that $\mathbf{U} \in B_{R_0}$. Then (4.10) and (4.11) yield the estimate

$$ \begin{aligned} |||\mathbf{h} + \mathbf{f}|||_{1,0,B_\varepsilon}^2 \le c \lambda^2 A^2 R_0^2 (1 + G_0^2 + G_0^4) (1 + R_0^2 + R_0^4) \\ + 2 |||\mathbf{f}|||_{1,0,B_\varepsilon}^2 \end{aligned} \quad (4.12) $$

<footnote n="3">³Without loss of generality we suppose that $\varepsilon \le 1$.</footnote>
and, by (3.15), we obtain

$$ \begin{aligned} & |||\mathcal{M}\mathbf{u}|||_{3,0,B_\epsilon}^2 + |||\nabla P|||_{1,0,B_\epsilon}^2 \\ & \leq c_1 \lambda^2 A^2 R_0^2 (1 + G_0^2 + G_0^4) (1 + R_0^2 + R_0^4) + c_2 |||\mathbf{f}|||_{1,0,B_\epsilon}^2. \end{aligned} \quad (4.13) $$

Put $M_0^2 = c_2 |||\mathbf{f}|||_{1,0,B_\epsilon}^2$ and $R_0^2 = 2M_0^2$. Suppose that

$$ \lambda^2 \le \frac{1}{2c_1 A^2 (1+G_0^2+G_0^4)(1+2c_2 |||\mathbf{f}|||_{1,0,B_\epsilon}^2 + 4c_2^2 |||\mathbf{f}|||_{1,0,B_\epsilon}^4)} = \lambda_*^2. \quad (4.14) $$

Then from (4.13) it follows that

$$ |||\mathcal{M}\mathbf{u}|||_{3,0,B_\epsilon}^2 \le R_0^2. $$

The last inequality implies that the operator $\mathcal{L}$ maps the closed bounded set $\mathcal{B}_{R_0} \subset W^{3,2}(B_\epsilon)$ into itself.

Let us show that $\mathcal{L}$ is a contraction in $H(B_\epsilon)$. Multiplying equations (4.9) by an arbitrary $\eta \in H(B_\epsilon)$ and integrating by parts we get

$$ \begin{aligned} \frac{\nu_0}{2} \int_{B_\epsilon} \nabla (\mathcal{M}\mathbf{U}) \cdot \nabla \eta \, dx = & -\lambda \int_{B_\epsilon} \nu(\dot{\gamma}(\mathbf{U} + \varepsilon\mathbf{g})) D(\mathbf{U} + \varepsilon\mathbf{g}) \cdot \nabla \eta \, dx \\ & +\lambda \int_{B_\epsilon} \nu(\dot{\gamma}(\varepsilon\mathbf{g})) D(\varepsilon\mathbf{g}) \cdot \nabla \eta \, dx + \int_{B_\epsilon} \mathbf{f} \cdot \eta \, dx. \end{aligned} \quad (4.15) $$

From (4.15) it follows that for any $\mathbf{U}_1, \mathbf{U}_2 \in \mathcal{B}_{R_0}$ the equality holds

$$ \begin{aligned} & \frac{\nu_0}{2} \int_{B_\epsilon} \nabla (\mathcal{M}\mathbf{U}_1 - \mathcal{M}\mathbf{U}_2) \cdot \nabla \eta \, dx \\ = & -\lambda \int_{B_\epsilon} \nu(\dot{\gamma}(\mathbf{U}_1 + \varepsilon\mathbf{g})) (D(\mathbf{U}_1) - D(\mathbf{U}_2)) \cdot \nabla \eta \, dx \\ & -\lambda \int_{B_\epsilon} (\nu(\dot{\gamma}(\mathbf{U}_1 + \varepsilon\mathbf{g})) - \nu(\dot{\gamma}(\mathbf{U}_2 + \varepsilon\mathbf{g}))) D(\mathbf{U}_2 + \varepsilon\mathbf{g}) \cdot \nabla \eta \, dx \\ = & J_1 + J_2. \end{aligned} \quad (4.16) $$

Since by (1.1),

$$ \begin{aligned} & |\nu(\dot{\gamma}(\mathbf{U}_1 + \varepsilon\mathbf{g})) - \nu(\dot{\gamma}(\mathbf{U}_2 + \varepsilon\mathbf{g}))|^2 \le \sup_{y} |\nabla_y \nu(y)|^2 |D(\mathbf{U}_1) - D(\mathbf{U}_2)|^2 \\ & \le A^2 |D(\mathbf{U}_1) - D(\mathbf{U}_2)|^2 \le cA^2 |\nabla \mathbf{U}_1 - \nabla \mathbf{U}_2|^2, \end{aligned} $$

we have

$$ |J_2| \le \frac{\nu_0}{8} \int_{B_\epsilon} |\nabla \eta|^2 dx + \frac{2c\lambda^2 A^2}{\nu_0} \int_{B_\epsilon} |\nabla \mathbf{U}_1 - \nabla \mathbf{U}_2|^2 |\nabla (\mathbf{U}_2 + \varepsilon \mathbf{g})|^2 dx $$

$$ \le \frac{\nu_0}{8} \int_{B_\epsilon} |\nabla \eta|^2 dx + \frac{2c\lambda^2 A^2}{\nu_0} \sup_{x \in B_\epsilon} |\nabla (\mathbf{U}_2 + \varepsilon \mathbf{g})|^2 \int_{B_\epsilon} |\nabla \mathbf{U}_1 - \nabla \mathbf{U}_2|^2 dx $$
$$ \leq \frac{\nu_0}{8} \int_{B_\epsilon} |\nabla \eta|^2 dx + \frac{2c\lambda^2 A^2}{\nu_0} (\epsilon^{2\alpha_1} |||\mathbf{U}_2|||_{3,0, B_\epsilon}^2 + G_0^2) \int_{B_\epsilon} |\nabla \mathbf{U}_1 - \nabla \mathbf{U}_2|^2 dx $$

$$ \leq \frac{\nu_0}{8} \int_{B_\epsilon} |\nabla \eta|^2 dx + \frac{2c\lambda^2 A^2}{\nu_0} (R_0^2 + G_0^2) \int_{B_\epsilon} |\nabla \mathbf{U}_1 - \nabla \mathbf{U}_2|^2 dx; $$

$$ |J_1| \leq \frac{\nu_0}{8} \int_{B_\epsilon} |\nabla \eta|^2 dx + \frac{2c\lambda^2 A^2}{\nu_0} \int_{B_\epsilon} |\nabla \mathbf{U}_1 - \nabla \mathbf{U}_2|^2 dx $$

Taking in (4.16) $\eta = \mathcal{M}\mathbf{U}_1 - \mathcal{M}\mathbf{U}_2$ we derive the inequality

$$ \frac{\nu_0}{2} \|\nabla(\mathcal{M}\mathbf{U}_1 - \mathcal{M}\mathbf{U}_2)\|_{L^2(B_\epsilon)}^2 \leq \frac{\nu_0}{4} \|\nabla(\mathcal{M}\mathbf{U}_1 - \mathcal{M}\mathbf{U}_2)\|_{L^2(B_\epsilon)}^2 $$

$$ + \lambda^2 \frac{c_3 A^2 [1 + R_0^2 + G_0^2]}{\nu_0} \|\nabla(\mathbf{U}_1 - \mathbf{U}_2)\|_{L^2(B_\epsilon)}^2. $$

Therefore,

$$ \|\nabla(\mathcal{M}\mathbf{U}_1 - \mathcal{M}\mathbf{U}_2)\|_{L^2(B_\epsilon)}^2 \leq \lambda^2 \frac{4c_3A^2[1 + R_0^2 + G_0^2]}{\nu_0^2} \|\nabla(\mathbf{U}_1 - \mathbf{U}_2)\|_{L^2(B_\epsilon)}^2. $$

Let

$$ \lambda_0^2 = \min \left\{ \lambda_*^2, \frac{\nu_0^2}{4c_3A^2[1 + R_0^2 + G_0^2]} \right\}. \qquad (4.17) $$

Then for any $\lambda \in (0, \lambda_0)$ the operator $\mathcal{M}$ is a contraction with the contraction factor

$$ q = \lambda^2 \frac{4c_3A^2[1 + R_0^2 + G_0^2]}{\nu_0^2} < 1 $$

and, by Theorem 3.10, there exists a unique fixed point $\mathbf{u}$ of the operator $\mathcal{M}$ which is a solution (together with the corresponding pressure function $p$) of problem (4.6). Estimate (4.7) for the fixed point $\mathbf{u}$ and the pressure $p$ follows from the fact that $\mathbf{u} \in B_R$ and inequality (4.13).

(ii) If the pressure $p$ satisfies the condition $\int_{B_\epsilon} p(x)dx = 0$, it can be represented in the form $p = \operatorname{div}\mathbf{w}$, where $\mathbf{w} \in \dot{W}^{1,2}(B_\epsilon)$ and

$$ \|\nabla \mathbf{w}\|_{L^2(B_\epsilon)} \leq c\epsilon^{-1} \|p\|_{L^2(B_\epsilon)}. $$
(see Lemma 3.6). Multiplying equations (4.6) by **w** and integrating by parts we derive

$$ \frac{\nu_0}{2} \int_{B_\epsilon} \nabla \mathbf{u} \cdot \nabla \mathbf{w} dx + \lambda \int_{B_\epsilon} \left( \nu(\dot{\gamma}(\mathbf{u} + \varepsilon \mathbf{g})) D(\mathbf{u} + \varepsilon \mathbf{g}) - \nu(\dot{\gamma}(\varepsilon \mathbf{g})) D(\varepsilon \mathbf{g}) \right) \cdot \nabla \mathbf{w} dx \\ - \int_{B_\epsilon} \mathbf{f} \cdot \mathbf{w} dx = \int_{B_\epsilon} p \mathrm{div} \mathbf{w} dx = \int_{B_\epsilon} |p|^2 dx. $$

Therefore,

$$ \int_{B_{\epsilon}} |p|^2 dx \le \|\mathbf{f}\|_{L^2(B_{\epsilon})} \|\mathbf{w}\|_{L^2(B_{\epsilon})} + \frac{\nu_0}{2} \|\nabla \mathbf{u}\|_{L^2(B_{\epsilon})} \|\nabla \mathbf{w}\|_{L^2(B_{\epsilon})} \\ + \lambda \left( \int_{B_{\epsilon}} |\nu(\dot{\gamma}(\mathbf{u} + \varepsilon \mathbf{g})) D(\mathbf{u} + \varepsilon \mathbf{g}) - \nu(\dot{\gamma}(\varepsilon \mathbf{g})) D(\varepsilon \mathbf{g})|^2 dx \right)^{1/2} \|\nabla \mathbf{w}\|_{L^2(B_{\epsilon})} \\ \le c\epsilon \|\mathbf{f}\|_{L^2(B_{\epsilon})} \|\nabla \mathbf{w}\|_{L^2(B_{\epsilon})} + \frac{\nu_0}{2} \|\nabla \mathbf{u}\|_{L^2(B_{\epsilon})} \|\nabla \mathbf{w}\|_{L^2(B_{\epsilon})} \\ + cA(1+G_0)^{1/2} \|\nabla \mathbf{u}\|_{L^2(B_{\epsilon})} \|\nabla \mathbf{w}\|_{L^2(B_{\epsilon})} \\ \le c \left( \|\mathbf{f}\|_{L^2(B_{\epsilon})} + \epsilon^{-1} \|\nabla \mathbf{u}\|_{L^2(B_{\epsilon})} \right) \|p\|_{L^2(B_{\epsilon})}, $$

and we obtain

$$ \|p\|_{L^2(B_\epsilon)}^2 \le c \left( \|\mathbf{f}\|_{L^2(B_\epsilon)}^2 + \epsilon^{-2} \|\nabla \mathbf{u}\|_{L^2(B_\epsilon)}^2 \right). $$

From the last inequality it follows that

$$ \|p\|_{L^2(B_\epsilon)}^2 \le c\epsilon^2 (\|\mathbf{f}\|_{1,0,B_\epsilon}^2 + \|\mathbf{u}\|_{3,0,B_\epsilon}^2) \le c\epsilon^2 \|\mathbf{f}\|_{1,0,B_\epsilon}^2. $$

□

# 5. NON-NEWTONIAN POISEUILLE FLOW

## 5.1. Existence of non-Newtonian Poiseuille flow with prescribed pressure slope

The non-Newtonian Poiseuille flow with the strain rate dependent viscosity was studied in the book [6] and recently in [27]. We will need below some extended versions of theorems proved there. Namely we will use the results on the regularity of this flow obtained in [26].

**Theorem 5.1.** (i) Let $\partial\sigma \in C^3$. For any $\alpha_0 > 0$ there exists $\lambda_0 = \lambda_0(\alpha_0)$ such that for all $\lambda \in (0, \lambda_0]$ and any $|\alpha| \le \alpha_0$ problem (1.5) admits a unique$^4$ solution $v_{P_\alpha} \in \dot{W}^{1,2}(\sigma) \cap W^{3,2}(\sigma)$. The solution $v_{P_\alpha}$ satisfies the estimate

$$ \|v_{P_\alpha}\|_{W^{3,2}(\sigma)} \le c|\alpha|, \quad \left| \int_{\sigma} v_{P_\alpha}(x')dx' \right| \le c|\alpha|, \qquad (5.1) $$

where the constant $c$ depends only on $\sigma, \nu_0, A$.

(ii) For any $F_0 > 0$ there exists $\lambda_1 = \lambda_1(F_0)$ such that for all $\lambda \in (0, \lambda_1]$ and every $F \in (-F_0, F_0)$ problem (1.5) admits a unique solution $(v_{P_\alpha}, \alpha)$ with $F_\sigma(\alpha) = \int_{\sigma} v_{P_\alpha}(x')dx' = F$. Moreover, the following estimates

$$ \|v_{P_\alpha}\|_{W^{3,2}(\sigma)} \le C|F|, \quad |\alpha| \le c|F|. \qquad (5.2) $$

<footnote>$^4$Here and below the uniqueness takes place only in some ball where the contraction principle is applied.</footnote>
hold.

(iii) Let $\alpha_1, \alpha_2$ be two real numbers such that $|\alpha_i| \le \alpha_0$, $i = 1, 2$, let $v_{P_{\alpha_1}}, v_{P_{\alpha_2}}$ be two solutions of problem (1.5). There exists $\lambda_3 = \lambda_3(\alpha_0)$ such that for all $\lambda \in (0, \lambda_3]$ the following estimates hold

$$ \left| \int_{\sigma} (v_{P_{\alpha_1}} - v_{P_{\alpha_2}}) dx' \right| \le c|\alpha_1 - \alpha_2|, \quad \|v_{P_{\alpha_1}} - v_{P_{\alpha_2}}\|_{W^{2,2}(\sigma)} \le c|\alpha_1 - \alpha_2|, \qquad (5.3) $$

where the constant $c$ depends only on $\sigma, \nu_0, A$.

The proof is given in [26]

## 5.2. Operator relating the pressure slope and the flux

Let us recall the notations from the Introduction: $F_\sigma(\alpha) = \int_\sigma v_{P_\alpha}(x')dx'$ is the flux corresponding to the pressure slope $-\alpha$, and $G(\alpha) = F_\sigma(\alpha) - \kappa\alpha$.

**Lemma 5.2.** For any $\alpha_0 > 0$ there exists a number $\lambda_2 = \lambda_2(\alpha_0) \le \lambda_1(\alpha_0)$ such that for all $\lambda \in (0, \lambda_2]$ the operator $\kappa^{-1}G(\alpha)$ is a contraction on the interval $[-\alpha_0, \alpha_0]$.

**Remark 5.3.** If the constant $\kappa^{-1}$ is replaced by another constant $K^{-1} > 0$ then for any $\alpha_0 > 0$ there exists a number $\lambda'_2 = \lambda'_2(\alpha_0) \le \lambda_1(\alpha_0)$ such that for all $\lambda \in (0, \lambda'_2]$ the operator $K^{-1}G(\alpha)$ is a contraction on the interval $[-\alpha_0, \alpha_0]$.

## 6. EQUATION ON THE GRAPH

Consider the following problem on the graph $\mathcal{B}$: given constants $F_l, l = N_1 + 1, \dots, N$, such that $\sum_{l=N_1+1}^N F_l = 0$ and constants $c_{lj}, l = 1, \dots, N_1$ (here for any $l$ subscript $j$ is such that the edges $e_j$ have an end point $O_l$), find a function $p$ which is affine with respect to $x_1^{(e_j)}$,

$$ p(x_1^{(e_j)}) = -s_j x_1^{(e_j)} + a_j, \qquad (6.1) $$

and such that

$$ \left\{ \begin{array}{l} -\frac{\partial}{\partial x_1^{(e_j)}} \left( F_{\sigma_j} \left( \frac{\partial p}{\partial x_1^{(e_j)}}(x_1^{(e_j)}) \right) \right) = 0, \quad x_1^{(e_j)} \in (0, |e_j|), \\ -\sum_{e_j: O_l \in e_j} F_{\sigma_j} \left( \frac{\partial p}{\partial x_1^{(e_j)}}(0) \right) = 0, \quad l = 1, \dots, N_1, \\ -F_{\sigma_j} \left( \frac{\partial p}{\partial x_1^{(e_j)}}(0) \right) = -F_l, \quad l = N_1 + 1, \dots, N, \quad O_l \in e_j, \\ p(x_1^{(e_j)} = 0) - p(x_1^{(e_s)} = 0) = c_{lj}, \quad e_j: O_l \in e_j, \quad l = 1, \dots, N_1, \\ p(O_N) = 0, \end{array} \right. \qquad (6.2) $$

where $e_s$ is a selected and fixed edge of the bundle.

This problem can be generalized as follows. Denote by $\mathcal{H}(\mathcal{B})$ the space of functions defined on the graph and belonging to $W^{1,2}(e_j)$ on every edge $e_j$ of the graph and vanishing at $O_N$. The norm in $\mathcal{H}(\mathcal{B})$ is defined by

$$ \|p\|_{\mathcal{H}(\mathcal{B})}^2 = \sum_{j=1}^{M} \|p\|_{W^{1,2}(e_j)}^2. $$
Given $F_l \in \mathbb{R}$, $l = 1, \dots, N$, $f^{(e_j)} \in L^2(e_j)$, $j = 1, \dots, M$, $c_{lj} \in \mathbb{R}$, $l = 1, \dots, N_1$, (for $j$ such that the edges $e_j$ have an end point $O_l$), and such that $\sum_{l=1}^N F_l = 0$ find a function $p \in \mathcal{H}(\mathcal{B})$ satisfying the equations

$$ egin{cases} 
-\frac{\partial}{\partial x_1^{(e_j)}} \left( F_{\sigma_j} \left( \frac{\partial p}{\partial x_1^{(e_j)}}(x_1^{(e_j)}) \right) \right) = f^{(e_j)}(x_1^{(e_j)}), \quad x_1^{(e_j)} \in (0, |e_j|), \\ 
-\sum_{e_j: O_l \in e_j} F_{\sigma_j} \left( \frac{\partial p}{\partial x_1^{(e_j)}}(0) \right) = -F_l, \quad l = 1, \dots, N_1, \\ 
-F_{\sigma_j} \left( \frac{\partial p}{\partial x_1^{(e_j)}}(0) \right) = -F_l, \quad l = N_1 + 1, \dots, N, \quad O_l \in e_j, \\ 
p(x_1^{(e_j)} = 0) - p(x_1^{(e_j)} = 0) = c_{lj}, \quad e_j: O_l \in e_j, \quad l = 1, \dots, N_1, \\ 
p(O_N) = 0. 
\end{cases} 
\quad (6.3) $$

in the sense of the following weak formulation:

$$ 
\begin{aligned} 
& \sum_{j=1}^{M} \int_{0}^{|e_j|} F_{\sigma_j} \left( \frac{\partial p}{\partial x_1^{(e_j)}}(x_1^{(e_j)}) \right) \frac{\partial q}{\partial x_1^{(e_j)}} dx_1^{(e_j)} + \sum_{l=1}^{N} F_l q(O_l) \\ 
& = \sum_{j=1}^{M} \int_{0}^{|e_j|} f^{(e_j)}(x_1^{(e_j)}) q(x_1^{(e_j)}) dx_1^{(e_j)} 
\end{aligned} 
\quad (6.4) $$

for all $q \in W^{1,2}(\mathcal{B})$ equal to 0 at $O_N$.

**Theorem 6.1.** There exists $\lambda_1$ such that $\forall \lambda \in [0, \lambda_1)$ problem (6.3) admits a unique weak solution from $\mathcal{H}(\mathcal{B})$. Let $f^{(e_j)(m)} \in L^2(e_j)$, $j = 1, \dots, M$, $F_l^{(m)}$, $c_{lj}^{(m)} \in \mathbb{R}$, $m = 1, 2$, be two sets of data and let $p^{(m)}$ be solutions of problem (6.3) corresponding to these data sets. Then there exists a constant $C$ depending on $\lambda, \sigma_j, \mathcal{B}$, such that

$$ 
\begin{aligned} 
\|p^{(1)} - p^{(2)}\|_{\mathcal{H}(\mathcal{B})}^2 \le C \left( \sum_{j=1}^{M} \|f^{(e_j)(1)} - f^{(e_j)(2)}\|_{L^2(e_j)}^2 \right. \\ 
\left. + \sum_{l=1}^{N} |F_l^{(1)} - F_l^{(2)}|^2 + \sum_{l=1}^{N_1} \sum_{e_j: O_l \in e_j} |c_{lj}^{(1)} - c_{lj}^{(2)}|^2 \right). 
\end{aligned} 
\quad (6.5) $$

Proof is given in [25].

**Corollary 6.2.** There exists $\lambda_1$ such that for all $\lambda \in [0, \lambda_1)$ problem (6.2) admits a unique solution from $\mathcal{H}(\mathcal{B})$. Let $F_l^{(m)} \in \mathbb{R}$, $l = N_1 + 1, \dots, N$, $c_{lj}^{(m)} \in \mathbb{R}$, $l = 1, \dots, N_1$, $m = 1, 2$, be two sets of data and let $p^{(m)}$ be solutions of problem (6.2) corresponding to these data sets. Then there exists a constant $C_p$ depending on $\lambda, \sigma_j, \mathcal{B}$, such that

$$ \|p^{(1)} - p^{(2)}\|_{\mathcal{H}(\mathcal{B})}^2 \le C \left( \sum_{l=N_1+1}^{N} |F_l^{(1)} - F_l^{(2)}|^2 + \sum_{l=1}^{N_1} \sum_{e_j: O_l \in e_j} |c_{lj}^{(1)} - c_{lj}^{(2)}|^2 \right). \quad (6.6) $$

## 7. SCALING OF THE NON-NEWTONIAN QUASI-POISEUILLE FLOW IN A THIN TUBE OF THE TUBE STRUCTURE. SCALING OF THE EQUATION ON THE GRAPH

In the introduction we have found the relations between the fluxes and Poiseuille velocities in the tubes with the section $\sigma$ and $\sigma_\varepsilon$:
$$v_{P, \frac{\alpha}{\varepsilon}}^{\varepsilon}(z') = \varepsilon v_{P\alpha}(\frac{z'}{\varepsilon}), \qquad (7.1)$$

$$F_{\sigma_{\varepsilon}}(\beta) = \varepsilon^{n} F_{\sigma}(\varepsilon\beta). \qquad (7.2)$$

Consider now the problem on the graph corresponding to the data of problem (1.2): the cross-sections $\sigma_{\varepsilon}^j$ and the given fluxes $F_l^{\varepsilon} = \int_{\gamma_{\varepsilon}^l} \varepsilon \mathbf{g}^l(\frac{x-O_l}{\varepsilon}) \cdot \mathbf{n}(x) dS = \varepsilon^n F^l,$

$$\left\{
\begin{array}{l}
-\frac{\partial}{\partial x_1^{(e_j)}} \left( F_{\sigma_{\varepsilon}^j} \left( \frac{\partial p^{\varepsilon}}{\partial x_1^{(e_j)}} (x_1^{(e_j)}) \right) \right) = 0, \quad x_1^{(e_j)} \in (0, |e_j|), \\
-\sum_{e_j: O_l \in e_j} F_{\sigma_{\varepsilon}^j} \left( \frac{\partial p^{\varepsilon}}{\partial x_1^{(e_j)}} (0) \right) = 0, \quad l = 1, \dots, N_1, \\
-F_{\sigma_{\varepsilon}^j} \left( \frac{\partial p^{\varepsilon}}{\partial x_1^{(e_j)}} (0) \right) = -F_l^{\varepsilon}, \quad l = N_1 + 1, \dots, N, \quad O_l \in e_j, \\
p^{\varepsilon}(x_1^{(e_j)} = 0) - p^{\varepsilon}(x_1^{(e_s)} = 0) = c_{lj}, \quad e_j: O_l \in e_j, \quad l = 1, \dots, N_1, \\
p^{\varepsilon}(O_N) = 0.
\end{array}
\right.
\qquad (7.3)$$

and it can be rescaled using the scalings of the pressure slope - flux relation:

$$\left\{
\begin{array}{l}
-\frac{\partial}{\partial x_1^{(e_j)}} \left( F_{\sigma_j}(\varepsilon \frac{\partial p^{\varepsilon}}{\partial x_1^{(e_j)}}(x_1^{(e_j)})) \right) = 0, \quad x_1^{(e_j)} \in (0, |e_j|), \\
-\sum_{e_j: O_l \in e_j} F_{\sigma_j}(\varepsilon \frac{\partial p^{\varepsilon}}{\partial x_1^{(e_j)}}(0)) = 0, \quad l = 1, \dots, N_1, \\
-F_{\sigma_j}(\varepsilon \frac{\partial p^{\varepsilon}}{\partial x_1^{(e_j)}}(0)) = -F^l, \quad l = N_1 + 1, \dots, N, \quad O_l \in e_j, \\
\varepsilon p^{\varepsilon}(x_1^{(e_j)} = 0) - \varepsilon p^{\varepsilon}(x_1^{(e_s)} = 0) = \varepsilon c_{lj}, \quad e_j: O_l \in e_j, \quad l = 1, \dots, N_1, \\
\varepsilon p^{\varepsilon}(O_N) = 0.
\end{array}
\right.
\qquad (7.4)$$

Notice that the operator $F_{\sigma_j}(\beta)$ is a nonlinear operator relating the pressure slope and the flux. Now making the change $p = \varepsilon p^{\varepsilon}$, we get a problem of (6.2) type where the left hand side does not depend on $\varepsilon$. If $p$ is an affine function, then

$$p^{\varepsilon}(x_1^{(e_j)}) = -s_j x_1^{(e_j)} / \varepsilon + a_j / \varepsilon. \qquad (7.5)$$

## 8. EXISTENCE, UNIQUENESS AND STABILIZATION OF A SOLUTION TO THE NON-NEWTONIAN FLOW EQUATIONS IN AN UNBOUNDED DOMAIN WITH CYLINDRICAL OUTLETS TO INFINITY

In this section we will recall the theorems from [26] which will be used in the construction of the boundary layer correctors for the asymptotic expansion of a solution of problem (1.2).
Consider in the domain $\Omega$ the steady state boundary value problem for the non-Newtonian fluid motion equations

$$ \left\{ \begin{array}{l} -\mathrm{div}((\nu_0 + \lambda\nu(\dot{\gamma}(\mathbf{v})))D(\mathbf{v})) + \nabla p = 0, \quad x \in \Omega, \\ \mathrm{div}\mathbf{v} = 0, \quad x \in \Omega, \\ \mathbf{v}|_{\partial\Omega} = 0. \end{array} \right. \qquad (8.1) $$

We look for the solution $\mathbf{v}$ having prescribed fluxes $F_j$ over the cross sections $\sigma_j$ of outlets to infinity:

$$ \int_{\sigma_j} \mathbf{v} \cdot \mathbf{n} dS = F_j, \quad j = 1, 2, \dots, J, \qquad (8.2) $$

where

$$ \sum_{j=1}^{J} F_j = 0. \qquad (8.3) $$

## 8.1. Existence and uniqueness of a solution

Consider the domain $\Omega \subset \mathbb{R}^n$ with $J$ cylindrical outlets to infinity. We assume that the boundary $\partial\Omega$ is $C^4$-regular. Consider in $\Omega$ the problem (8.1), (8.2), (8.3). Denote $\mathbb{F} = \sqrt{\sum_{j=1}^{J} F_j^2}$. By Theorem 5.1, for any set of fluxes $(F_1, \dots, F_J)$ such that $\mathbb{F} \le F_0$, there is a number $\lambda_{00}$ depending on $F_0$, such that for every $\lambda \in (0, \lambda_{00})$ there exist $J$ pressure slopes $\alpha_j$ and corresponding $J$ quasi-Poiseuille flows $\mathbf{V}_{P_{\alpha_j}}(x) = (v_{P_{\alpha_j}}(x'), 0, \dots, 0)^T \in W^{3,2}(\sigma_j)$, defined in cylinders $\{x^{(j)} \in \mathbb{R}^n, x^{(j)'} \in \sigma_j, x_1^{(j)} \in \mathbb{R}\}, j = 1, \dots, J$, such that $F_{\sigma_j}(\alpha_j) = F_j$.

We define the cut-off functions $\chi_j$ associated to each outlet $\Omega^j$ as $C^4$-smooth functions vanishing everywhere in $\Omega$ except for the outlet $\Omega^j$, where they depend on the local longitudinal variable $x_1^{(j)}$ only, are equal to zero if $x_1^{(j)} < 1$ and equal to one if $x_1^{(k)} > 2$. Put

$$ \mathbf{V}_{\chi} = \sum_{j=1}^{J} \chi_j \mathbf{V}_{P_{\alpha_j}}, \quad P_{\chi} = -\sum_{j=1}^{J} \chi_j \alpha_j x_1^{(e_j)}. $$

It is easy to see that for $h$ given by the formula

$$ h(x) = \mathrm{div} \, \mathbf{V}_{\chi}(x) = \sum_{j=1}^{J} \chi_j'(x_1^{(j)}) v_{P_{\alpha_j}}(x^{(j)'}) $$
$$ \mathrm{supp} \, h \subset \overline{\Omega^{(2)} \setminus \Omega^{(1)}}. \qquad (8.4) $$

Moreover, from the condition $\sum_{j=1}^{J} F_j = 0$ it follows that

$$ \int_{\Omega^{(2)}} h(x) dx = 0. $$
Finally, estimates (5.1) and (5.2) yield

$$ \|h\|_{W^{2,2}(\Omega^{(3)})} \le c \sum_{j=1}^{J} \|v_{P_{\alpha_j}}\|_{W^{2,2}(\sigma_j)} \le c \mathbb{F} \quad (8.5) $$

Since $h \in W^{2,2}(\Omega^{(3)})$, by results in [8], there exists a vector field $\mathbf{W} \in \dot{W}^{1,2}(\Omega^{(3)}) \cap W^{3,2}(\Omega^{(3)})$ such that

$$ \mathrm{div} \, \mathbf{W}(x) = h(x), $$

and

$$ \|\mathbf{W}\|_{W^{3,2}(\Omega^{(3)})} \le c \|h\|_{W^{2,2}(\Omega^{(3)})} \le c \mathbb{F}. \quad (8.6) $$

Moreover, since $\mathrm{supp} h \subset \bar{\Omega}^{(2)}$, $\mathbf{W}$ can be constructed such that

$$ \mathrm{supp} \, \mathbf{W} \subset \bar{\Omega}^{(3)}. \quad (8.7) $$

Extend the functions $\mathbf{W}$ and $\mathbf{V}_\chi$ by zero into the whole $\Omega$ and set

$$ \hat{\mathbf{V}}_{\chi}(x) = \mathbf{W}(x) + \mathbf{V}_{\chi}(x). \quad (8.8) $$

Then,

$$ \mathrm{div} \, \hat{\mathbf{V}}(x) = 0, \quad \hat{\mathbf{V}}(x)|_{\partial \Omega} = 0, \quad \int_{\sigma_j} \hat{\mathbf{V}}(x) \cdot \mathbf{n}(x) \, ds = F_j, \quad j = 1, \dots, J, $$

and for $x \in \Omega^j \setminus \Omega^{(3)}$ the vector-field $\hat{\mathbf{V}}(x)$ coincides with the velocity part $\mathbf{V}_{P_{\alpha_j}}(x^{(j)})$ of the corresponding Poiseuille flow.

By denoting in (1.2)

$$ \mathbf{v} = \mathbf{u} + \hat{\mathbf{V}}_{\chi}, \quad p = q + \mathcal{P}_{\chi}, \quad (8.9) $$

where $\mathcal{P}_\chi = \sum_{j=1}^J \chi_j \alpha_j x_1^j$, we obtain the following problem

$$ \left\{ \begin{array}{l} -\mathrm{div} \left[ (\nu_0 + \lambda \nu(\dot{\gamma}(\mathbf{u} + \hat{\mathbf{V}}_{\chi}))D(\mathbf{u} + \hat{\mathbf{V}}_{\chi})) + \nabla(q + \mathcal{P}_{\chi}) \right] = \mathbf{f} \quad \text{in } \Omega, \\ \mathrm{div} \mathbf{u} = 0 \quad \text{in } \Omega, \\ \mathbf{u} = 0 \quad \text{on } \partial\Omega, \\ \int_{\sigma_j} \mathbf{u} \cdot \mathbf{n} \, dS = 0, \quad j = 1, \dots, J. \end{array} \right. \quad (8.10) $$

**Theorem 8.1.** For any $f_0 > 0$ and $\mathbb{F}_0 > 0$ there exist numbers $\Lambda_0 = \Lambda_0(F_0, f_0) > 0$ and $\beta_* > 0$ such that $\forall \lambda \in (0, \Lambda_0]$, $\forall \beta \in (0, \beta_*]$ and for any $\mathbf{f} \in W_{\beta}^{1,2}(\Omega)$ satisfying $\|\mathbf{f}\|_{W_{\beta}^{1,2}(\Omega)} \le f_0$ and any set $(F_1, \dots, F_J)$ with
$\mathbb{F}^2 = \sum_{j=1}^{J} F_j^2 \le F_0^2$ the problem (8.1), (8.2), (8.3) has a unique solution $(\mathbf{v}, p)^5$ admitting the representation (8.9) with $\mathbf{u} \in W_{\beta}^{3,2}(\Omega)$, $\nabla q \in W_{\beta}^{1,2}(\Omega)$, $\int_{\Omega^{(3)}} q(x)dx = 0$. The following estimate

$$ \|\mathbf{u}\|_{W_{\beta}^{3,2}(\Omega)}^2 + \|\nabla q\|_{W_{\beta}^{1,2}(\Omega)}^2 \le c \left( \|\mathbf{f}\|_{W_{\beta}^{1,2}(\Omega)}^2 + \mathbb{F}^2 \right) \quad (8.11) $$

holds. Moreover, there exist constants $q_1, q_2, \dots, q_J$ such that

$$ \begin{aligned} & \int_{\Omega_0} |q(x)|^2 dx + \sum_{j=1}^{J} \int_{\Omega_j} \exp(2\beta x_1^{(j)}) |q(x) - q_j|^2 dx \\ & \le c \int_{\Omega} E_{\beta}(x) |\nabla q(x)|^2 dx \le c \left( \|\mathbf{f}\|_{W_{\beta}^{1,2}(\Omega)}^2 + \mathbb{F}^2 \right). \end{aligned} \quad (8.12) $$

The proof is given in [26].

Below we will say that the function $q \in L_{\text{loc}}^2(\Omega)$ exponentially stabilizes to constants $q_1, q_2, \dots, q_J$ at infinity if

$$ \int_{\Omega^j} \exp(2\beta x_1^{(j)}) |q(x) - q_j|^2 dx < \infty, \quad j = 1, 2, \dots, J, $$

for some $\beta > 0$. The space of such functions we denote $\tilde{L}_{\text{loc}}^2(\Omega)$.

## 8.2. Continuity of the solution with respect to data of the problem

Assume that we have two sets of fluxes $(F_1^{(1)}, \dots, F_1^{(J)})$ and $(F_1^{(2)}, \dots, F_1^{(J)})$ satisfying the condition (8.3), and two functions $\mathbf{f}^{(1)}, \mathbf{f}^{(2)} \in W_{\beta}^{1,2}(\Omega)$. Let $\hat{\mathbf{V}}_\chi^{(1)}$ and $\hat{\mathbf{V}}_\chi^{(2)}$ be flux carriers corresponding to fluxes $(F_1^{(1)}, \dots, F_1^{(J)})$ and $(F_1^{(2)}, \dots, F_1^{(J)})$, respectively (see formula (8.8)). Denote by $(\mathbf{u}^{(1)}, q^{(1)})$ and $(\mathbf{u}^{(2)}, q^{(2)})$ the solutions of problem (8.10) corresponding to the flux carriers $\hat{\mathbf{V}}_\chi^{(1)}, \hat{\mathbf{V}}_\chi^{(2)}$ and the right hand sides $\mathbf{f}^{(1)}, \mathbf{f}^{(2)}$. Assume that

$$ \mathbb{F}^{(i)} \le F_0, \quad \|\mathbf{f}^{(i)}\|_{W_{\beta}^{1,2}(\Omega)} \le f_0, \quad i = 1, 2. \quad (8.13) $$

Denote

$$ Q = \sum_{j=1}^{J} |F_j^{(1)} - F_j^{(2)}|^2 + \|\mathbf{f}^{(1)} - \mathbf{f}^{(2)}\|_{W_{\beta}^{1,2}(\Omega)}^2. $$

**Theorem 8.2.** There exists $\Lambda = \Lambda_1(F_0, f_0)$ and $\beta_*$ such that for $\forall \lambda \in (0, \Lambda_1]$, $\forall \beta \in (0, \beta_*]$ and sufficiently small $Q$ for arbitrary $\mathbf{f}^{(i)}$ and $(F_1^{(i)}, \dots, F_J^{(i)})$, $i = 1, 2$, satisfying (8.13), the following estimate holds:

$$ \|\mathbf{u}^{(1)} - \mathbf{u}^{(2)}\|_{W^{2,2}(\Omega)}^2 + \|\nabla(q^{(1)} - q^{(2)})\|_{L^2(\Omega)}^2 \le cQ |\ln Q|^2. \quad (8.14) $$

<footnote>$^5$The uniqueness takes place only in some ball where the contraction principle is applied and we have in mind the uniqueness only for solutions admitting the representation (8.9). Moreover, as usual, the pressure $p$ is unique up to an additive constant.</footnote>
Moreover, there exist constants $(q_1^{(1)}, \dots, q_J^{(1)})$ and $(q_1^{(2)}, \dots, q_J^{(2)})$ such that if $q^{(1)}$ and $q^{(2)}$ are normalized by the conditions $\int_{\Omega^{(3)}} q^{(1)}(x)dx = \int_{\Omega^{(3)}} q^{(2)}(x)dx = 0$, then

$$ \sum_{j=1}^{J} \int_{\Omega_j} \exp(2\beta x_1^{(j)}) |q^{(m)}(x) - q_j^{(m)}|^2 dx \le c(F_0 + f_0), \quad m = 1, 2, \qquad (8.15) $$

and

$$ \sum_{j=1}^{J} |q_j^{(1)} - q_j^{(2)}|^2 \le cQ |\ln Q|^2. \qquad (8.16) $$

The proof is given in [26].

## 9. CONSTRUCTION OF AN ASYMPTOTIC APPROXIMATION OF THE SOLUTION

The algorithm of the construction of an asymptotic approximation is a sequence of iterations where at each step we solve a problem on the graph and a set of boundary layer problems (see [21]). Namely, the asymptotic expansion consists of two addends. The first addend is a quasi-Poiseuille flow multiplied by a cut-off function vanishing within some neighbourhood of the nodes/vertices $O_l$. Substituting this first addend in problem (1.2), we derive N boundary layer problems which are obtained by scaling of problem (8.10) in the neighbourhood of nodes and vertices. These two addends are matched via the constants $c_{lj}$ in the problem on the graph (7.3) and the constants $\tilde{c}_{lj}$ to which stabilizes at infinity the pressure in (8.10). Constants $c_{lj}$ and $\tilde{c}_{lj}$ should be the same. By an iterative process this condition will be achieved with the accuracy $O(\varepsilon^J)$.

**Step 0.1.** Given $F_l$, $l = N_1 + 1, \dots, N$, we first solve on the graph the problem of order zero

$$ \left\{ \begin{array}{l} -\frac{\partial}{\partial x_1^{(e_j)}} \left( F_{\sigma_j} \left( \frac{\partial p_0}{\partial x_1^{(e_j)}} (x_1^{(e_j)}) \right) \right) = 0, \quad x_1^{(e_j)} \in (0, |e_j|), \\ -\sum_{e_j: O_l \in e_j} F_{\sigma_j} \left( \frac{\partial p_0}{\partial x_1^{(e_j)}} (0) \right) = 0, \quad l = 1, \dots, N_1, \\ -F_{\sigma_j} \left( \frac{\partial p_0}{\partial x_1^{(e_j)}} (0) \right) = -F^l, \quad l = N_1 + 1, \dots, N, \quad O_l \in e_j, \\ p_0(x_1^{(e_j)} = 0) - p_0(x_1^{(e_s)} = 0) = 0, \quad e_j: O_l \in e_j, \quad l = 1, \dots, N_1, \\ p_0(O_N) = 0. \end{array} \right. \qquad (9.1) $$

and define $p_0^\varepsilon = \varepsilon^{-1}p_0$,

$$ p_0^\varepsilon(x_1^{(e_j)}) = -s_{j,0}x_1^{(e_j)}/\varepsilon + a_{j,0}/\varepsilon, \quad v_{P,\alpha_{j,0}^\varepsilon}^\varepsilon(x^{(e_j)}) = \varepsilon v_{P,s_{j,0}}(x^{(e_j)})/\varepsilon. \qquad (9.2) $$

Multiplying by a cut-off function we introduce for $x_1^{(e_j)} \in (0, |e_j|/2)$ the functions

$$ \begin{aligned} \mathbf{v}_0^a &= v_{P, \alpha_{j,0}^\varepsilon}^\varepsilon(x^{(e_j)}) \zeta\left(\frac{x_1^{(e_j)}}{3r\varepsilon}\right) \mathbf{e}_j; \\ p_0^a &= (p_0^\varepsilon(x_1^{(e_j)}) - p_0^\varepsilon(x_1^{(e_s)} = 0)) \zeta\left(\frac{x_1^{(e_j)}}{3r\varepsilon}\right) + p_0^\varepsilon(x_1^{(e_s)} = 0), \end{aligned} \qquad (9.3) $$

where $\zeta(t) = 0$ if $|t| \le 1$, $\zeta(t) = 1$ if $|t| \ge 2$, and $\zeta \in C^3(\mathbb{R})$ and $\mathbf{e}_j$ is the direction vector of the edge $e_j$.
**Step 0.2.** Denote $\bar{x} = x - O_l$, $\mathbf{V}_{\zeta}^{\epsilon}(\bar{x}) = \sum_{j:O_l \in e_j} \zeta(\frac{x_1^{(e_j)}}{3r\epsilon})v_{0,P,\alpha_{j,0}^{\epsilon}}(x^{(e_j)\prime})\mathbf{e}_j$, $P_{\zeta}^{\epsilon}(\bar{x}) = \sum_{j:O_l \in e_j} \zeta(\frac{x_1^{(e_j)}}{3r\epsilon})s_j x_1^{(e_j)}/\epsilon$. Substituting (9.3) into (1.2) and formally extending the branches to infinity, we get for each $O_l$ the following problem posed in $\Omega_l^{\epsilon}$: find $\mathbf{N}_0^{\epsilon} \in \dot{W}^{1,2}(\Omega_l^{\epsilon})$ and $P_0^{\epsilon} \in \tilde{L}_{loc}^2(\Omega_l^{\epsilon})$, such that

$$ 
\left\{ 
\begin{aligned}
& -\operatorname{div}((\nu_0 + \lambda\nu(\dot{\gamma}(\mathbf{N}_0^{\epsilon}(\bar{x}) + \mathbf{V}_{\zeta}^{\epsilon}(\bar{x})))D(\mathbf{N}^{\epsilon}(\bar{x}) + \mathbf{V}_{\zeta}^{\epsilon}(\bar{x}))) \\
& \qquad + \nabla(P_0^{\epsilon} - P_{\zeta}^{\epsilon}(\bar{x})) = 0, \quad \bar{x} \in \Omega_l^{\epsilon}, \\
& \operatorname{div}(\mathbf{N}_0^{\epsilon}(\bar{x}) + \mathbf{V}_{\zeta}^{\epsilon}(\bar{x})) = 0, \quad \bar{x} \in \Omega_l^{\epsilon}, \\
& \mathbf{N}_0^{\epsilon} = 0, \quad \bar{x} \in \partial\Omega_l^{\epsilon},
\end{aligned}
\right. 
\quad (9.4) $$

and so in dilated variables $\xi = \frac{x-O_l}{\epsilon}$ in the unbounded domain $\Omega_l$: find $\mathbf{N}_0 \in \dot{W}^{1,2}(\Omega_l)$ and $P_0 \in \tilde{L}_{loc}^2(\Omega_l)$, such that

$$ 
\left\{ 
\begin{aligned}
& -\operatorname{div}_{\xi}((\nu_0 + \lambda\nu(\dot{\gamma}_{\xi}(\mathbf{N}_0(\xi) + \mathbf{V}_{\zeta}^0(\xi))))D_{\xi}(\mathbf{N}_0(\xi) + \mathbf{V}_{\zeta}^0(\xi))) \\
& \qquad + \nabla_{\xi}(P_0(\xi) - P_{\zeta}^0(\xi)) = 0, \quad \xi \in \Omega_l, \\
& \operatorname{div}_{\xi}(\mathbf{N}_0(\xi) + \mathbf{V}_{\zeta}^0(\xi)) = 0, \quad \xi \in \Omega_l, \\
& \mathbf{N}_0(\xi) = 0, \quad \xi \in \partial\Omega_l,
\end{aligned}
\right. 
\quad (9.5) $$

where

$$ \mathbf{V}_{\zeta}^{0}(\xi) = \sum_{j:O_l \in e_j} \zeta(\frac{\xi_1^{(e_j)}}{3r})v_{P,s_{j,0}}(\xi^{(e_j)\prime})\mathbf{e}_j, P_{\zeta}^{0}(\xi) = \sum_{j:O_l \in e_j} \zeta(\frac{\xi_1^{(e_j)}}{3r})s_{j,0}\xi_1^{(e_j)}, $$

$\operatorname{div}_{\xi}, \dot{\gamma}_{\xi}, D_{\xi}, \nabla_{\xi}$ are operators written for $\xi$ variable. Notice that this problem is independent of $\epsilon$, $\mathbf{N}_0^{\epsilon}(x - O_l) = \epsilon\mathbf{N}_0(\frac{x-O_l}{\epsilon})$, $P_0^{\epsilon}(x - O_l) = P_0(\frac{x-O_l}{\epsilon})$. Let $P_0$ tends to zero at the outlet $\Pi_{\infty}^{(e_s)}$ corresponding to the selected edge $e_s$ of the bundle $B_l$, and denote by $\tilde{c}_{lj}^0$ the constants which are the limits of $P_0$ in the outlets $\Pi_{\infty}^{(e_j)}$, $j \neq s$. Here and below for vertices $O_l$ condition $\mathbf{N}_k = 0$ on $\gamma^l$ is replaced by $\mathbf{N}_k = \mathbf{g}^l$; in this case there is only one outlet and the pressure tends there to zero. Then we pass to the order one.

**Step 1.1.** Given $F_l$, $l = N_1 + 1, \dots, N$ and constants $\tilde{c}_{lj}^0$ (all constants are independent of $\epsilon$) we solve again problem on the graph of the first order, that is

$$ 
\left\{ 
\begin{array}{l}
-\frac{\partial}{\partial x_1^{(e_j)}}\left(F_{\sigma_j}\left(\frac{\partial p_1}{\partial x_1^{(e_j)}}(x_1^{(e_j)})\right)\right) = 0, \quad x_1^{(e_j)} \in (0, |e_j|), \\
-\sum_{e_j: O_l \in e_j} F_{\sigma_j}\left(\frac{\partial p_1}{\partial x_1^{(e_j)}}(0)\right) = 0, \quad l = 1, \dots, N_1, \\
-F_{\sigma_j}\left(\frac{\partial p_1}{\partial x_1^{(e_j)}}(0)\right) = -F^l, \quad l = N_1 + 1, \dots, N, \quad O_l \in e_j, \\
p_1(x_1^{(e_j)} = 0) - p_1(x_1^{(e_s)} = 0) = \epsilon \tilde{c}_{lj}^0, \quad e_j: O_l \in e_j, \quad l = 1, \dots, N_1, \\
p_1(O_N) = 0.
\end{array}
\right. 
\quad (9.6) $$

and define $p_1^{\epsilon} = \epsilon^{-1}p_1$,

$$ p_1^{\epsilon}(x_1^{(e_j)}) = -s_{j,1}x_1^{(e_j)}/\epsilon + a_{j,1}/\epsilon, \quad (9.7) $$
and $v_{P,\alpha_{j,1}}^{\varepsilon}(x^{(e_j)}) = \varepsilon v_{P,s_{j,1}}(x^{(e_j)}/\varepsilon)$. Consider the functions

$$
\begin{align}
\mathbf{v}_1^a &= V_{P,\alpha_{j,1}^{\varepsilon}}^{\varepsilon}(x^{(e_j)}) \zeta\left(\frac{x_1^{(e_j)}}{3r\varepsilon}\right) \mathbf{e}_j; \\
p_1^a &= \left(p_1^{\varepsilon}(x_1^{(e_j)}) - p_1^{\varepsilon}(x_1^{(e_s)} = 0)\right) \zeta\left(\frac{x_1^{(e_j)}}{3r\varepsilon}\right) + p_1^{\varepsilon}(x_1^{(e_s)} = 0), \tag{9.8} \\
x_1^{(e_j)} &\in (0, |e_j|/2).
\end{align}
$$

Applying the estimates of Theorem 5.1, Corollary 6.2 and Theorem 8.2, we get the estimates for the differences
$p_1 - p_0$, $s_{j,1} - s_{j,0}$, $v_{P,s_{j,1}}^{\varepsilon} - v_{P,s_{j,0}}^{\varepsilon}$:

$$
\begin{align}
\|p_1 - p_0\|_{\mathcal{H}(\mathcal{B})}^2 &\le C \sum_{l=1}^{N_1} \sum_{e_j: O_l \in e_j} |\varepsilon(c_{lj}^1 - c_{lj}^0)|^2, \nonumber \\
\sum_{j=1}^M |s_{j,1} - s_{j,0}|^2 &\le C \sum_{l=1}^{N_1} \sum_{e_j: O_l \in e_j} |\varepsilon(c_{lj}^1 - c_{lj}^0)|^2, \tag{9.9} \\
\sum_{j=1}^M \|v_{P,s_{j,1}} - v_{P,s_{j,0}}\|_{H^1(\sigma_j)}^2 &\le C \sum_{l=1}^{N_1} \sum_{e_j: O_l \in e_j} |\varepsilon(c_{lj}^1 - c_{lj}^0)|^2, \nonumber
\end{align}
$$

where $C$ is a constant independent of $\varepsilon$ and $c_{lj}^0 = 0$, $c_{lj}^1 = \tilde{c}_{lj}^0$.

**Step 1.2.** Solve the following problem in dilated variables $\xi = \frac{x-O_l}{\varepsilon}$ in the unbounded domain $\Omega_l$: find $\mathbf{N}_1 \in \mathbf{H}^1(\Omega_l)$ and $P_1 \in \tilde{\mathcal{L}}_{\text{loc}}^2(\Omega_l)$, such that

$$
\left\{
\begin{array}{l}
-\operatorname{div}_{\xi}((\nu_0 + \lambda\nu(\dot{\gamma}_{\xi}(\mathbf{N}_1(\xi) + \mathbf{V}_{\zeta}^1(\xi))))D_{\xi}(\mathbf{N}_1(\xi) + \mathbf{V}_{\zeta}^1(\xi))) \\
\quad + \nabla_{\xi}(P_1(\xi) - P_{\zeta}^1(\xi)) = 0, \quad \xi \in \Omega_l, \\
\operatorname{div}_{\xi}(\mathbf{N}_1(\xi) + \mathbf{V}_{\zeta}^1(\xi)) = 0, \quad \xi \in \Omega_l, \\
\mathbf{N}_1(\xi) = 0, \quad \xi \in \partial\Omega_l.
\end{array}
\right.
\tag{9.10}
$$

Here and below

$$
\mathbf{V}_{\zeta}^{1}(\xi) = \sum_{j:O_l \in e_j} \zeta\left(\frac{\xi_{1}^{(e_j)}}{3r}\right) v_{0,P,s_{j,1}}(\xi^{(e_j)}) \mathbf{e}_j
$$

and

$$
P_{\zeta}^{1}(\xi) = \sum_{j:O_l \in e_j} \zeta\left(\frac{\xi_1^{(e_j)}}{3r}\right) s_{j,1} \xi_1^{(e_j)}, \quad k = 1, \dots, J.
$$

As in the previous step we assume that $P_1$ tend to zero in the selected outlet $\Pi_{\infty}^{(e_s)}$, and by $\tilde{c}_{lj}^1$ are denoted the constants which are the limits of $P_1$ at the outlets $\Pi_{\infty}^{(e_j)}$.

Applying the estimates of Theorem 8.2, we get that

$$
|\tilde{c}_{lj}^{1} - \tilde{c}_{lj}^{0}| \leq C\varepsilon |\ln \varepsilon|^2.
$$

Assume that we have constructed functions $p_k$, constants $s_{j,k}$, $a_{j,k}$, functions $v_{P,s_{j,k}}$ and for all $l = 1, \dots, N$,
functions $\mathbf{N}_k$, $P_k$ and constants $c_{lj}^k = \tilde{c}_{lj}^{k-1}$ such that

$$
|c_{lj}^k - c_{lj}^{k-1}| \leq C \varepsilon^{k-1} |\ln \varepsilon|^{2(k-1)}
$$
and

$$
\begin{align}
\|p_k - p_{k-1}\|_{\mathcal{H}(B)}^2 &\le C \sum_{l=1}^{N_1} \sum_{e_j:O_l \in e_j} |\varepsilon(c_{lj}^k - c_{lj}^{k-1})|^2, \nonumber \\
\sum_{j=1}^{M} |s_{j,k} - s_{j,k-1}|^2 &\le C \sum_{l=1}^{N_1} \sum_{e_j:O_l \in e_j} |\varepsilon(c_{lj}^k - c_{lj}^{k-1})|^2, \tag{9.11} \\
\sum_{j=1}^{M} \|v_{P,s_{j,k}} - v_{P,s_{j,k-1}}\|_{W^{1,2}(\sigma^j)}^2 &\le C \sum_{l=1}^{N_1} \sum_{e_j:O_l \in e_j} |\varepsilon(c_{lj}^k - c_{lj}^{k-1})|^2, \nonumber
\end{align}
$$

where $\sqrt{\sum_{l=1}^{N_1} \sum_{e_j:O_l \in e_j} |\varepsilon(c_{lj}^k - c_{lj}^{k-1})|^2} = O(\varepsilon^k |\ln \varepsilon|^{2(k-1)})$. Define

$$
Q_k = \sum_{l=1}^{N_1} \sum_{e_j:O_l \in e_j} |\varepsilon(c_{lj}^k - c_{lj}^{k-1})|^2 = O(\varepsilon^{2k} |\ln \varepsilon|^{4(k-1)}).
$$

Then

$$
\begin{align}
\|\mathbf{N}_k - \mathbf{N}_{k-1}\|_{W^{1,2}(\Omega)}^2 &\le C Q_k |\ln Q_k| = O(\varepsilon^{2k} |\ln \varepsilon|^{4k-3}), \nonumber \\
\|P_k - P_{k-1} - \sum_{j:O_l \in e_j} \zeta(\frac{\xi_{1j}^{(e_j)}}{3r})(\tilde{c}_{lj}^k - \tilde{c}_{lj}^{k-1})\|_{L^2(\Omega)}^2 &\le C Q_k |\ln Q_k| \nonumber \\
&= O(\varepsilon^{2k} |\ln \varepsilon|^{4k-3}), \tag{9.12} \\
|\tilde{c}_{lj}^k - \tilde{c}_{lj}^{k-1}|^2 &\le C Q_k |\ln Q_k|^2 = O(\varepsilon^{2k} |\ln \varepsilon|^{4k-2}). \nonumber
\end{align}
$$

Define now $c_{lj}^{k+1} = \tilde{c}_{lj}^k$. So,

$$
|c_{lj}^{k+1} - c_{lj}^{k}| = O(\varepsilon^k |\ln \varepsilon|^{2k-1}) = O(\varepsilon^k |\ln \varepsilon|^{2k}).
$$

**Step k+1.1.** Given $F_l$, $l = N_1 + 1, \dots, N$ and constants $\tilde{c}_{lj}^k$ (these constants are uniformly bounded with respect to $\varepsilon$), we solve the problem on the graph of the $(k+1)$ order, that is

$$
\left\{
\begin{array}{l}
-\displaystyle\frac{\partial}{\partial x_1^{(e_j)}}\left(F_{\sigma_j}\left(\frac{\partial p_{k+1}}{\partial x_1^{(e_j)}}(x_1^{(e_j)})\right)\right) = 0, \quad x_1^{(e_j)} \in (0, |e_j|), \\
-\displaystyle\sum_{e_j:O_l \in e_j} F_{\sigma_j}\left(\frac{\partial p_{k+1}}{\partial x_1^{(e_j)}}(0)\right) = 0, \quad l = 1, \dots, N_1, \\
-F_{\sigma_j}\left(\frac{\partial p_{k+1}}{\partial x_1^{(e_j)}}(0)\right) = -F^l, \quad l = N_1 + 1, \dots, N, \quad O_l \in e_j, \\
p_{k+1}(x_1^{(e_j)} = 0) - p_{k+1}(x_1^{(e_s)} = 0) = \varepsilon \tilde{c}_{lj}^k, \quad O_l \in e_j, \quad l = 1, \dots, N_1, \\
p_{k+1}(O_N) = 0.
\end{array}
\right.
\tag{9.13}
$$

and define $p_{k+1}^{\varepsilon} = \varepsilon^{-1} p_{k+1}$,

$$
\begin{align*}
p_{k+1}^{\varepsilon}(x_1^{(e_j)}) &= -s_{j,k+1}x_1^{(e_j)}/\varepsilon + a_{j,k+1}/\varepsilon, \\
v_{0,P,\alpha_{j,k+1}^{\varepsilon}}^{\varepsilon}(x^{(e_j)\prime}) &= \varepsilon v_{0,P,s_{j,k+1}}(x^{(e_j)\prime}/\varepsilon).
\end{align*}
$$
Applying Theorem 5.1, Theorem 6.1 and Theorem 8.2, we get the estimates for the differences $p_{k+1} - p_k$, $s_{j,k+1} - s_{j,k}$, $v_{P,s_{j,k+1}}^{\epsilon} - v_{P,s_{j,k}}^{\epsilon}$:

$$ 
\begin{aligned}
\|p_{k+1} - p_k\|_{\mathcal{H}(\mathcal{B})}^2 &\le C \sum_{l=1}^{N_1} \sum_{e_j: O_l \in e_j} |\epsilon(\tilde{c}_{lj}^k - \tilde{c}_{lj}^{k-1})|^2 \\
&\le C \epsilon^2 Q_k |\ln Q_k| = O(\epsilon^{2(k+1)} |\ln \epsilon|^{4k}), \\
\sum_{j=1}^{M} |s_{j,k+1} - s_{j,k}|^2 &= O(\epsilon^{2(k+1)} |\ln \epsilon|^{4k}), \\
\sum_{j=1}^{M} \|v_{P,s_{j,k+1}} - v_{P,s_{j,k}}\|_{W^{1,2}(\sigma^j)}^2 &= O(\epsilon^{2(k+1)} |\ln \epsilon|^{4k}).
\end{aligned} 
\quad (9.14) $$

Define now $c_{lj}^{k+1} = \tilde{c}_{lj}^k$.

**Step k+1.2.** Solve problem in dilated variables $\xi = \frac{x-O_l}{\epsilon}$ in the unbounded domain $\Omega_l$: find $N_{k+1} \in W^{1,2}(\Omega_l)$ and $P_{k+1} \in \tilde{L}_{loc}^2(\Omega_l)$, such that

$$ 
\left\{ 
\begin{array}{l}
-\mathrm{div}_{\xi}((\nu_0 + \lambda\nu(\dot{\gamma}_{\xi}(\mathbf{N}_{k+1}(\xi) + \mathbf{V}_{\zeta}^{k+1}(\xi))))D_{\xi}(\mathbf{N}_{k+1}(\xi) + \mathbf{V}_{\zeta}^{k+1}(\xi))) \\
+\nabla_{\xi}(P_{k+1}(\xi) - P_{\zeta}^{k+1}(\xi)) = 0, \quad \xi \in \Omega_l, \\
\mathrm{div}_{\xi}(\mathbf{N}_{k+1}(\xi) + \mathbf{V}_{\zeta}^{k+1}(\xi)) = 0, \quad \xi \in \Omega_l, \\
\mathbf{N}_{k+1}(\xi) = 0, \quad \xi \in \partial\Omega_l.
\end{array}
\right. 
\quad (9.15) $$

Let $\tilde{c}_{lj}^{k+1}$ be the constants which are the limits of $P_{k+1}(\xi)$ at the outlets $\Pi_{\infty}^{(e_j)}$ for the edges of the bundle $\mathcal{B}_l$. Applying Theorem 8.2, we get the estimates

$$ 
\begin{aligned}
\|\mathbf{N}_{k+1} - \mathbf{N}_k\|_{W^{1,2}(\Omega)}^2 &\le C Q_{k+1} |\ln Q_{k+1}| = O(\epsilon^{2k+2} |\ln \epsilon|^{4k+1}), \\
\|P_{k+1} - P_k - \sum_{j:O_l \in e_j} \zeta(\frac{\xi^{(e_j)}}{3r}) (\tilde{c}_{lj}^{k+1} - \tilde{c}_{lj}^k)\|_{L^2(\Omega)}^2 &\le C Q_{k+1} |\ln Q_{k+1}| = O(\epsilon^{2k+2} |\ln \epsilon|^{4k+1}), \\
|\tilde{c}_{lj}^{k+1} - \tilde{c}_{lj}^k|^2 &\le C Q_{k+1} |\ln Q_{k+1}|^2 = O(\epsilon^{2k+2} |\ln \epsilon|^{4k+2}).
\end{aligned} 
\quad (9.16) $$

So, the algorithm is based on the fact that if the difference $c_{lj}^k - c_{lj}^{k-1}$ is of order $\epsilon^{k-1} |\ln \epsilon|^{2(k-1)}$, then $\|p_k - p_{k-1}\|_{\mathcal{H}(\mathcal{B})}$ is of order $\epsilon^k |\ln \epsilon|^{2(k-1)}$ due to the scaling between $p_k^{\epsilon}$ and $p_k = \epsilon p_k^{\epsilon}$; then $\|v_{p,s_{j,k}} - v_{p,s_{j,k-1}}\|_{W^{1,2}(\sigma^j)}$ is of the same order as the pressure, then $\tilde{c}_{lj}^k - \tilde{c}_{lj}^{k-1}$ is of order $\epsilon^k |\ln \epsilon|^{2k}$. So, $c_{lj}^{k+1} - c_{lj}^k = O(\epsilon^k |\ln \epsilon|^{2k})$.

Let us define the asymptotic expansion of the solution in each part of the domain $\mathcal{B}_{\epsilon}$ corresponding to the bundle $B_{O_l}$ truncated at the distance $|e_j|/2$ for every edge of the bundle $e_j$:

$$ 
\begin{aligned}
\mathbf{v}_J^a &= v_{P,\alpha_{j,j}^{\epsilon}}^{\epsilon}(x^{(e_j)}) \zeta(\frac{x_1^{(e_j)}}{3r\epsilon}) \mathbf{e}_j \mathcal{X}_j + \epsilon \mathbf{N}_J(\frac{x-O_l}{\epsilon}) (1 - \zeta(\frac{6x_1^{(e_j)}}{|e_j|}) \mathcal{X}_j) \\
&\quad + \Phi_{e_j}(x) \mathcal{X}_j; \\
p_J^a &= (p_J^{\epsilon}(x_1^{(e_j)}) - p_J^{\epsilon}(x_1^{(e_j)} = 0)) \zeta(\frac{x_1^{(e_j)}}{3r\epsilon}) \mathcal{X}_j + p_J^{\epsilon}(x_1^{(e_j)} = 0) \\
&\quad + (P_J(\frac{x-O_l}{\epsilon}) - \tilde{c}_{lj}^J \zeta(\frac{x_1^{(e_j)}}{3r\epsilon}) \mathcal{X}_j) (1 - \zeta(\frac{6x_1^{(e_j)}}{|e_j|}) \mathcal{X}_j) \\
&\quad + (\tilde{c}_{lj}^J - c_{lj}^J) \zeta(\frac{x_1^{(e_j)}}{3r\epsilon}) (1 - \zeta(\frac{6x_1^{(e_j)}}{|e_j|})) \mathcal{X}_j, \quad x \in B_{O_l}, x_1^{(e_j)} \in (0, |e_j|/2),
\end{aligned} 
\quad (9.17) $$

where $\mathcal{X}_j$ is the characteristic function of $\Pi_{\epsilon}^{(e_j)}$, $\Phi_{e_j}$ is a vector valued function with the support in $\{x^{(e_j)} : x_1^{(e_j)} \in (\frac{1}{7}|e_j|, \frac{2}{5}|e_j|), x^{\prime(e_j)} \in \sigma^j\}$, vanishing at the lateral boundary of the cylinder $\Pi_{\epsilon}^{(e_j)}$ and such that
within the domain $\hat{\Upsilon}_{\varepsilon}^{(e_j)} = \Pi_{\varepsilon}^{(e_j)} \cap \{x^{(e_j)} : x_1^{(e_j)} \in (\frac{1}{8}|e_j|, \frac{1}{2}|e_j|), x'(e_j) \in \sigma^j\}$,

$$ \mathrm{div}\Phi_{e_j} = -\mathrm{div}\{\varepsilon\mathbf{N}_J(\frac{x-O_l}{\varepsilon})(1-\zeta(\frac{6x_1^{(e_j)}}{|e_j|}))\}\mathcal{X}_j $$

(see Lem. 3.7). The factor $1 - \zeta(\frac{6x_1^{(e_j)}}{|e_j|})$ is introduced to remove the tails of $\mathbf{N}_J$ far from the nodes. Notice that near the note the divergence is equal zero due to the construction of $\mathbf{N}_J$ (see the second equation in (9.15)).

Note that $\mathbf{v}_J^a$ satisfies exactly the boundary condition and $\mathrm{div}\mathbf{v}_J^a = 0$. Let us calculate the discrepancy in the equations (1.2). The couple $(v_{P,\alpha_{j,J}}^{\varepsilon}(x^{(e_j)})\mathbf{e}_j, p_J^{\varepsilon}(x_1^{(e_j)}))$ is an exact solution of the equations within the cylinders $\Pi_{\varepsilon}^{(e_j)}$, the couple

$$ (v_{P,\alpha_{j,J}}^{\varepsilon}(x^{(e_j)})'\zeta(\frac{x_1^{(e_j)}}{3r\varepsilon})\mathbf{e}_j + \varepsilon\mathbf{N}_J(\frac{x-O_l}{\varepsilon}), $$

$$ (p_J^{\varepsilon}(x_1^{(e_j)}) - p_J^{\varepsilon}(x_1^{(e_j)} = 0))\zeta(\frac{x_1^{(e_j)}}{3r\varepsilon}) + p_J^{\varepsilon}(x_1^{(e_j)} = 0) + P_J(\frac{x-O_l}{\varepsilon}) $$

is an exact solution of the equations within the bundle $B_{O_l}$. This pressure term can be rewritten as

$$ (p_J^{\varepsilon}(x_1^{(e_j)}) - p_J^{\varepsilon}(x_1^{(e_j)} = 0) - c_{lj}^J)\zeta(\frac{x_1^{(e_j)}}{3r\varepsilon}) + p_J^{\varepsilon}(x_1^{(e_j)} = 0) + (P_J(\frac{x-O_l}{\varepsilon}) - \tilde{c}_{lj}^J\zeta(\frac{x_1^{(e_j)}}{3r\varepsilon})) \\ + \tilde{c}_{lj}^J\zeta(\frac{x_1^{(e_j)}}{3r\varepsilon}), $$

or,

$$ (p_J^{\varepsilon}(x_1^{(e_j)}) - p_J^{\varepsilon}(x_1^{(e_j)} = 0))\zeta(\frac{x_1^{(e_j)}}{3r\varepsilon}) + p_J^{\varepsilon}(x_1^{(e_j)} = 0) + (P_J(\frac{x-O_l}{\varepsilon}) - \tilde{c}_{lj}^J\zeta(\frac{x_1^{(e_j)}}{3r\varepsilon})) \\ + (\tilde{c}_{lj}^J - c_{lj}^J)\zeta(\frac{x_1^{(e_j)}}{3r\varepsilon}). $$

So, the residual $\mathbf{r}_{\varepsilon}(x)$ is generated by the difference between constants $\tilde{c}_{lj}^J$ and $c_{lj}^J$ (the last term in the approximation of the pressure), by the cut-off factor $1 - \zeta(\frac{6x_1^{(e_j)}}{|e_j|})\mathcal{X}_j$ in the boundary layer term and by the divergence corrector $\Phi_{e_j}$ (because the multiplication of the boundary layer by the cut off function slightly violates the "divergence free" equation):

$$ \mathbf{r}_{\varepsilon} = -\mathrm{div}((\nu_0 + \lambda\nu(\dot{\gamma}(\Psi)))D\Psi + \nabla\psi), \qquad (9.18) $$
where

$$ 
\begin{aligned}
\Psi(x) &= \varepsilon N_J\left(\frac{x-O_l}{\varepsilon}\right)\left(1-\zeta\left(\frac{6x_1^{(e_j)}}{|e_j|}\right)\right)\chi_j + \Phi_{e_j}(x)\chi_j; \\
\psi(x) &= \left(P_J\left(\frac{x-O_l}{\varepsilon}\right) - \tilde{c}_{lj}^J\zeta\left(\frac{x_1^{(e_j)}}{3r\varepsilon}\right)\chi_j\right)\left(1-\zeta\left(\frac{6x_1^{(e_j)}}{|e_j|}\right)\right)\chi_j \\
&\quad +\left(\tilde{c}_{lj}^J - c_{lj}^J\right)\zeta\left(\frac{x_1^{(e_j)}}{3r\varepsilon}\right)\left(1-\zeta\left(\frac{6x_1^{(e_j)}}{|e_j|}\right)\right)\chi_j.
\end{aligned} 
\quad (9.19) $$

The difference between constants $\tilde{c}_{lj}^J$ and $c_{lj}^J$ is of order $O(\varepsilon^J |\ln \varepsilon|^{2J+2})$. The multiplication of the term

$$ P_J\left(\frac{x - O_l}{\varepsilon}\right) - \tilde{c}_{lj}^J \zeta\left(\frac{x_1^{(e_j)}}{3r\varepsilon}\right) $$

by the cut-off factor $1 - \zeta(\frac{6x_1^{(e_j)}}{|e_j|})\chi_j$ brings a residual of order $O(e^{-c/\varepsilon})$ in $W^{1,2}$ norm with some positive $c$ independent of $\varepsilon$, due to the exponential stabilization of $P_J(\xi)$ to the constant $\tilde{c}_{lj}^J$ as $\xi$ tends to infinity in an outlet (so it is exponentially small in subdomain where $\zeta$ changes its value from one to zero; only in this subdomain the boundary layer correctors do not satisfy exactly the equations). Let us evaluate

$$ |\mathrm{div}((\nu_0 + \lambda\nu(\dot{\gamma}(\Psi)))D\Psi| \le c|\nabla^2\Psi| + c\lambda A|\nabla^2\Psi|(|\nabla\Psi| + 1), $$

$$ |\nabla \mathrm{div}((\nu_0 + \lambda\nu(\dot{\gamma}(\Psi)))D\Psi| \le c(|\nabla^3\Psi| + c\lambda A|\nabla^2\Psi|^2)(|\nabla\Psi| + 1), $$

and so,

$$ \begin{aligned} \|\mathbf{r}_{\varepsilon}\|_{H^1(B_{\varepsilon})} &\le c(\|(|\nabla^3\Psi| + |\nabla^2\Psi| + |\nabla^2\Psi|^2)(1 + |\nabla\Psi|)\|_{L^2(B_{\varepsilon})} \\ &\le c(\varepsilon^{-(n-4)/2}|||\Psi|||_{3,0,B_{\varepsilon}} + 1)^2|||\Psi|||_{3,0,B_{\varepsilon}} = O(e^{-c/\varepsilon}); \end{aligned} \quad (9.20) $$

$$ |||\nabla\psi|||_{1,0,B_{\varepsilon}} = O(\varepsilon^{J-2+(n-1)/2} |\ln \varepsilon|^{2J+2}). $$

Taking together these estimates we conclude that in $W^{1,2}(B_{\varepsilon})$ norm the residual has order $O(\varepsilon^{J-2} |\ln \varepsilon|^{2J+2})\varepsilon^{(n-1)/2}$.

Applying Theorem 4.1, we prove the following theorem.

**Theorem 9.1.** Assume that $g^j \in W^{5/2,2}(\gamma^j)$ satisfy condition (1.3) and let $\nu$ satisfy conditions (1.1). There exists $\hat{\lambda} > 0$ independent of small parameter such that for all $\lambda \in [0, \hat{\lambda})$ there exists a solution of problem (1.2), such that

$$ \begin{aligned} \|\mathbf{v}_J^a - \mathbf{v}_\varepsilon\|_{W^{1,2}(B_\varepsilon)} &\le C\varepsilon^{J-2}|\ln \varepsilon|^{2J+2}, \\ \|p_J^a - p_\varepsilon\|_{L^2(B_\varepsilon)} &\le C\varepsilon^{J-3}|\ln \varepsilon|^{2J+2} \end{aligned} \quad (9.21) $$

with a constant $C$ independent of $\varepsilon$.

**Remark 9.2.** The obtained estimates can be applied to justify the method of asymptotic partial decomposition of the domain for the non-Newtonian flows in thin tube structures. This method was tested in [17] however the question of the reconstruction of the pressure is still an open question in the case of the non-Newtonian flows.
## 10. CONCLUSION

The paper considers the stationary non-Newtonian Stokes system of equations in a thin tube structure modeling a network of blood vessels with the no slip boundary condition. The viscosity depends on the shear rate. The small parameter is the ratio of the radius of a vessel to its length. The complete asymptotic expansion of the solution is constructed and justified by the error estimates. These estimates evaluate the limitations of the asymptotic approximation, however, the numerical experiments often confirm wider frames for the asymptotic theory. The leading term of the asymptotic expansion is presented by (1.15) corresponding to the Poiseuille-like flow in the inner part of the vessels corrected with some boundary layer functions in the neighborhood of the junction zones. The Poiseuille part of the leading term is defined by a non-linear elliptic problem on the graph with the Kirchhoff type junction conditions at the nodes for the macroscopic pressure (1.12). Comparing this approximation to the one-dimensional models of the blood flow [5], we notice that the model in [5] is essentially non-stationary, having hyperbolic type; it is suitable for the description of the flows with high Reynolds' number when the inertial term dominates over the viscous term, while (1.15) is more advantageous for the quasi-stationary flows with dominating viscous term and modest Reynolds' number. Also the derivation of (1.15) starts with the non-Newtonian version of the Stokes system of equations while [5] derives the one-dimensional model directly from the conservation laws.

*Acknowledgements.* This project has received funding from European Social Fund (project No 09.3.3-LMT-K-712-17-0003) under grant agreement with the Research Council of Lithuania (LMTLT).

*Conflict of interest.* The authors declare that they have no conflict of interest.

## REFERENCES

[1] R.A. Adams, Sobolev Spaces. Academic Press, New York (1975).

[2] S. Agmon, A. Douglis and L. Nirenberg, Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions. I. *Commun. Pure Appl. Math.* **12** (1959) 623-727; II, *Commun. Pure Appl. Math.* **17** (1964) 35-92.

[3] M.E. Bogovskii, Solutions of some problems of vector analysis related to operators div and grad, *Proc. Semin. S.L. Sobolev* **1** (1980) 5-40 (in Russian).

[4] R. Bounoiu, A. Gaudiello and A. Leopardi, Asymptotic analysis of a Bingham fluid in a thin T-like shaped structure. *J. Math Pures Appl.* **123** (2019) 148-166.

[5] L. Formaggia, D. Lamponi and A. Quarteroni, One dimensional models for blood flow in arteries. *J. Eng. Math.* **47** (2003) 251-276.

[6] P. Galdi, R. Rannacher, A.M. Robertson and S. Turek, Hemodynamical Flows, Modeling, Analysis and Simulation. Oberwolfach Seminars, V.37, Birkhäuser, Basel, Boston, Berlin (2008).

[7] P. Galdi, An Introduction to the Mathematical Theory of Navier-Stokes Equations. Springer (1994).

[8] L.V. Kapitanskii and K. Pileckas, Certain problems of vector analysis, Zapiski Nauchn. Sem. LOMI, 1984, 138, 65-85 Engl. transl. *Sov. Math.* **32** (1986) 469-483.

[9] N. Kloviene and K. Pileckas, Non-stationary Poiseuille-type solutions for the second-grade fluid flow. *Lithuanian Math. J.* **52** (2012) 155-171.

[10] N. Kloviene and K. Pileckas, The second grade fluid flow problem in an infinite pipe. *Asymptotic Anal.* **83** (2013) 237-262.

[11] O.A. Ladyzhenskaya, The Mathematical Theory of Viscous Incompressible Fluid. Gordon and Breach (1969).

[12] O.A. Ladyzhenskaya and V.A. Solonnikov, On some problems of vector analysis and generalized formulations of boundary value problems for the Navier-Stokes equations. *Zapiski Nauchn. Sem. LOMI* **59** (1976) 81-116. English Transl.: *J. Sov. Math.* **10** (1978) 257-285.

[13] O.A. Ladyzhenskaya and N.N. Ural'ceva, Linear and Quasi-Linear Elliptic Equations. Academic Press, New York (1968).

[14] E.M. Landis, Second Order Elliptic and Parabolic Equations. Nauka, Moscow (1971).

[15] E. Marušić-Paloka, Steady flow of a non-Newtonian fluid in unbounded channels and pipes. *Math. Models Methods Appl. Sci.* **10** (2000) 1425-1445.

[16] E. Marušić-Paloka and I. Pažanin, A note on Kirchhoff junction rule for power-law fluids. *Zeitschrift für Naturforschung A* **70** (2015) 695-702.

[17] A. Nachit, G. Panasenko and A.M. Zine, Asymptotic partial domain decomposition in thin tube structures: numerical experiments. *Int. J. Multiscale Comput. Eng.* **11** (2013) 407-441.

[18] S.A. Nazarov and B.A. Plamenevskiy, Elliptic Problems in Domains with Piecewise Smooth Boundary. Nauka, Moscow (1991).

[19] O.A. Oleinik and G.A. Yosifian, On the asymptotic behaviour at infinity of solutions in linear elasticity. *Arch. Rat. Mech. Anal.* **78** (1982) 29-53.
[20] G. Panasenko, Asymptotic expansion of the solution of Navier-Stokes equation in a tube structure, *C.R. Acad. Sci. Paris 326*, Série IIb (1998) 867-872.

[21] G. Panasenko, Multi-scale Modeling for Structures and Composites. Springer, Dordrecht (2005).

[22] G. Panasenko and K. Pileckas, Divergence equation in thin-tube structure. *Applicable Anal. 94* (2015) 1450-1459.

[23] G. Panasenko and K. Pileckas, Flows in a tube structure: equation on the graph. *J. Math. Phys. 55* (2014) 081505.

[24] G. Panasenko and K. Pileckas, Asymptotic analysis of the non-steady Navier-Stokes equations in a tube structure. I. The case without boundary layer-in-time. *Nonlinear Anal. Series A, Theory Methods Appl. 122* (2015) 125-168.

[25] G. Panasenko, K. Pileckas and B. Vernescu, Steady state non-Newtonian flow in thin tube structure: equation on the graph. *Algebra Anal. 33* (2021) 197-214.

[26] G. Panasenko, K. Pileckas and B. Vernescu, Steady state non-Newtonian flow with strain rate dependent viscosity in domains with cylindrical outlets to infinity. *Nonlinear Anal.: Modeling Controle 26* (2021) 1166-1199.

[27] G. Panasenko and B. Vernescu, Non-Newtonian flows in domains with non-compact boundaries. *Nonlinear Anal. A 183* (2019) 214-229.

[28] K. Pileckas, Navier-Stokes system in domains with cylindrical outlets to infinity. Leray's problem, *Handbook of Mathematical Fluid Dynamics*, 4. Elsevier (2007), Chap. 8, p. 445-647.

[29] K. Pileckas, Weighted $L^q$-solvability of the steady Stokes system in domains with incompact boundaries. *Math. Models Methods Appl. Sci. 6* (1996) 97-136.

[30] K. Pileckas, On the non-stationary linearized Navier-Stokes problem in domains with cylindrical outlets to infinity. *Math. Ann. 332* (2005) 395-419.

[31] K. Pileckas, A. Sequeira and J.H. Videman, Steady flows of viscoelastic fluids in domains with outlets to infinity. *J. Math. Fluid Mech. 2* (2000) 185-218.

[32] K. Rajagopal and A. Gupta, On a class of exact solutions to the equations of motion of a 2D grade fluids. *J. Eng. Sci. 19* (1981) 1009-1014.

[33] K. Rajagopal, A note on unsteady unidirectional flows of a non-Newtonian fluid. *Int. J. Nonlinear Mech. 17* (1982) 369-373.

Subscribe to Open (S2O)
A fair and sustainable open access model

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

This journal is currently published in open access under a **Subscribe-to-Open** model (S2O). S2O is a transformative model that aims to move subscription journals to open access. Open access is the free, immediate, online availability of research articles combined with the rights to use these articles fully in the digital environment. We are thankful to our subscribers and sponsors for making it possible to publish this journal in open access, free of charge for authors.

**Please help to maintain this journal in open access!**

Check that your library subscribes to the journal, or make a personal donation to the S2O programme, by contacting
**subscribers@edpsciences.org**

More information, including a list of sponsors and a financial transparency report, available at:
**https://www.edpsciences.org/en/maths-s2o-programme**