Numerical Study of the Equation on the Graph for the Steady State non-Newtonian Flow in Thin Tube Structure

. The dimension reduction for the viscous flows in thin tube structures leads to equations on the graph for the macroscopic pressure with Kirchhoff type junction conditions in the vertices. Non-Newtonian rheology of the flow generates nonlinear equations on the graph. A new numerical method for second order nonlinear differential equations on the graph is introduced and numerically tested.


Introduction and main definitions
The Newtonian flows in a network of tubes (tube structures) were considered in [2,5,13,15,16,18,19].These domains are finite connected unions of thin cylinders (in the two-dimensional case, respectively, thin rectangles).Each tube structure is represented by its graph: if the diameter of cross-sections of cylinders tends to zero, then the tubes degenerate and "tend" to the edges of the graph.In the case of the stationary Navier-Stokes equations, set in a thin structure, for the pressure one can derive second order ordinary differential equations at the edges of the graph and some Kirchhoff-type junction conditions at the nodes (see [2,13,15,17]).The Kirchhoff junction conditions were introduced for some physiological applications in classical works [10,11].
In the non-stationary case for the initial-boundary value problem for the Navier-Stokes system, the equation on the graph becomes non-stationary and non-local in time (see [19]).For the periodic in time setting, it was studied in [20].In [7], a linear stationary equation in an infinite fractal type graph is considered.
Nonlinear equations on the graph are generated by the problems in thin tube structures for viscous flows with non-Newtonian rheology.The case of the Stokes equations with the shear rate dependent viscosity was considered in the set of three papers [21,22,23], where the complete asymptotic expansion of the solution was constructed (also see [3,4] for the Bingham law rheology and [9] for the power law rheology).In particular, in [22] a nonlinear equation on the graph is studied.Its unknown function is the macroscopic pressure in the network of thin pipes.This problem is crucial for the construction of asymptotic approximations of the solution of the Stokes type equation in tube structures because it gives the leading term.These problems on the graph need fast effective solvers.The present paper introduces a new numerical method of its solution, formulates the convergence theorem and provides numerical tests.
Other one-dimension models for the blood flow applications were derived from the general conservation laws in [6,12,25].

Thin tube structure and its graph
For the reader's convenience we remind the definitions of the tube structure and its graph given in [13].3, and e 1 , e 2 , ..., 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 ∈ {1, ..., N }, i j ̸ = k j ).All points O i are supposed to be the ends of some segments e j .The segments e j are called edges of the graph.A point O i is called a node, if it is the common end of at least two edges, and O i is called a vertex, if it is the end of the only one edge.Any two edges e j and e i can intersect only at the common node.The set of vertices is supposed to be non-empty.Denote B = ∪ M j=1 e j the union of edges and assume that B is a connected set.The graph G is defined as the collection of nodes, vertices and edges (see Figure 1).
Let e be some edge, e = O i O j .Consider two Cartesian coordinate systems in R n .The first one has the origin in O i and the axis O i x  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 σ j ⊂ R n−1 containing the origin O and having C 4 -smooth boundary ∂σ j , j = 1, ..., M .For every edge e j = e and associated σ j = σ (e) we denote by σ (e) ε the set x (e)′ ∈ R n−1 : x (e)′ /ε ∈ σ (e) , and by Π (e) ε the cylinder ("tube") where n ), |e| is the length of the edge e and ε > 0 is a small parameter.Notice that the edges e j and Cartesian coordinates of nodes and vertices O j , as well as the domains σ j , do not depend on ε.
Definition 2. By a tube structure we call the following domain (see Figure 2) Suppose that it is a connected set and that the boundary The fourth order smoothness of the boundary is required to prove the existence of the solution to the non-Newtonian problem, see [23].Consider a node or a vertex O l and all edges e j having O l as one of their end points.We call a bundle of edges B l the union of all these edges, i.e., B l = ∪ j:O l ∈ej e j .
3 Motivation: a non-Newtonian flow in network of tubes Let n = 2, 3, ν 0 , λ > 0 be positive constants.Let ν be a bounded where C is a positive constant.
In the tube structure B ε , consider the steady state boundary value problem for the non-Newtonian fluid motion equations where D(v) is the strain rate matrix with the elements Assume that the fluid velocity g at the boundary ∂B ε has the following structure: g = 0 everywhere on ∂B ε except for the set γ N1+1 ε , ..., γ N ε , where where Let e = e Oj be the edge with the end at the vertex O j and let x (e) be the Cartesian coordinates corresponding to the origin O j and the edge e, i.e., x (e) = P (e) (x − O j ), P (e) is the orthogonal matrix relating the global coordinates x with the local ones x (e) .Denote g (e) = P (e) g j .Let where n is the unit outward (with respect to B ε ) normal vector to γ j ε , y (e) = x (e) /ε, ĝj (y (e) ) = g j ((P (e) ) T y (e) ) (here "T" denotes the transposition), F j does not depend on ε.Assume the compatibility condition for the flow rates F j : J j=1 F j = 0.An asymptotic analysis of this problem with small parameter ε leads to a nonlinear elliptic problem on the graph [21,22,23].The solution of this problem determines the leading term of the pressure in problem (3.1).
Let us introduce this problem on the graph.

Equation on the graph
For any edge e with associated cross-section σ (e) where v Pα is a solution of the boundary value problem , and α is a given pressure slope.In what follows, we omit the subscript ε.
With the above notations, consider the following problem on the graph B: given constants F l , l = N 1 + 1, ..., N , such that N l=N1+1 F l = 0 and constants c lj , l = 1, ..., 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 The fourth condition means that the solution p may have prescribed jumps c lj at the nodes O l ∈ e j , l = 1, ..., N 1 , so that for a fixed node O l at the ends of different edges e j (when the local variable x (ej ) 1 is equal to 0) the solution p has its value which differs from the value of p at the end of the selected edge e s and the difference is equal to c lj .For the leading term c lj = 0, and so p is continuous on the graph.However for the further terms c lj may be different from zero.
Remark 1.In [21,22,23] it is proved that for any α 0 > 0, C > 0 there exists λ 0 > 0 such that for |λ| ≤ λ 0 , CG σ (e) is a contraction on the interval [−α 0 , α 0 ], and there exists a solution to (4.2), unique in some ball B R = {q ∈ H O N (B), ∥q∥ H(B) ≤ R}.Here H(B) is the space of functions defined on the graph, belonging to H 1 (e j ) on every edge e j of the graph, and H O N (B) is the subspace of functions from H(B) vanishing at O N ; the space H(B) is supplied with the norm defined by the relation: We will also consider the subspaces H 1 (B) of H(B) and consisting of continuous on B functions, without gaps at the nodes.
The proof of the existence of the solution to problem (4.2) is based on the fixed point method applied in some ball B R in the space H(B).Every iteration is related to the following linear problem.
For f = 0) − p(x Function p is continuous on B and satisfies the following variational formulation: )dx Here It is supposed that the following compatibility condition holds: Treating the left-hand side of (4.4) as an inner product in the space H 1 O N (B) and its right-hand side as a bounded linear functional we apply the Riesz theorem and prove that this problem admits a unique solution.Moreover, taking the test function q = p and applying the Cauchy-Buniakovsky-Schwarz inequality, we prove that there exists a constant C B depending on B only, such that For more details, see [17].
Remark 2. Note that if for any edge e, f Lp(x This problem is an elementary iteration to solve problem (4.2).Applying Remark 1, we get that there exists λ 0 such that for all |λ| ≤ λ 0 , L is a contraction in H L (B) ∩ B R .Therefore, problem (4.2) can be solved by the fixed point iterations converging for any initial approximation from B R , for example, p 0 = 0.The iterations are defined as follows: (5.1)Each piecewise-affine function p ∈ H L (B) ∩ B R is in one-to-one correspondence with the set of values of p in the ends of the edges e, so that (5.1) is equivalent to a linear system of equations with inversible matrix M. So, it can be solved by applying the LU -factorization of M (calculated only once for the whole chain of iterations) and then solving the two standard triangular systems of equations at each step of iterations.However, for large N and sparse matrix M, the GMRES method for each iteration may be more effective.

Fixed point method with preconditioner
Consider a more robust iterative method: given a positive number τ (further called "artificial time"), calculate recurrently ∈ (0, |e j |), (5.2) Let us analyse for what values of τ this method converges.Let λ be a positive constant such that C B G σ (e) is a contraction with factor ρ. Let p be the exact solution of problem (4.2).Then, evidently, the difference q k = p k − p is a solution of the problem consider a real life blood rheology corresponding to the non-Newtonian Carreau law [8] ν 0 + λν( γ) = 0.00345 + 0.05255(1 + (3.313 γ) 2 ) −0.3216 Pa • s.
To define the function F σ (e) ε (α), we assume that the cylinders of the tube structure have a round shape of cross-sections, so that σ  In Figure 3, one can see that the level surfaces of the velocity are parallel everywhere except for the small zones near the bases of the cylinder (boundary layer zones).The smallness of the boundary layer zones confirms the asymptotic analysis developed in [23] in the case when the thin tube structure is presented by one thin cylinder.According to this analysis on the cross-section in the middle of the tube, we obtain the solution of the corresponding problem (4.1). Figure 4 shows the computed relation between the flux and the pressure slope for several values of the pressure slope.
To fit it for a tube σ ε with diameter equal to 0.01ε m, the following property (see [22]) is employed:  The diameters of the vessels are assumed to be equal to 10 −3 m and their lengths to 10 −2 m.For realistic fluxes, the following boundary conditions are used: 1 below shows values of nodes/vertices O i computed with the iterative algorithm.An exact solution is presented as well, which was found using the MATLAB nonlinear equations solver Fsolve.Note that these results are subject to one predetermined unknown (O 1 in this case), which is required for a unique solution to exist.

A realistic graph example
Next, the iterative algorithm was tested on a realistic graph of blood vessels.A segment of network from [24] was used with realistic lengths, while the diameters were simplified to all equal 10 −4 m. Figure 6 below depicts the graph of this vessel network.

1
has the direction of the ray [O i O j ); the second one has the origin in O j and the opposite direction, i.e.O ix(e) 1 is directed over the ray [O j O i ).

(4. 3 )
Let us describe a change of the unknown function p which allows to reduce problem (4.3) to a problem of the same type but for an unknown function continuous on the graph.Let ζ be a C 3 -smooth function defined on [0, +∞) equal to zero on [0, 1/6] and equal to one on [1/3, +∞).Introducing on every edge e j function 1−ζ x (ej ) 1 /|e j | c lj we change the unknown function p on the edge e j (j ̸ = s) by p = p − 1 − ζ x (ej ) 1 /|e j | c lj .We get for the new unknown function p(e) .Equation (4.3) with new right-hand sides

1 =
const, then on each edge e the solution p is an affine function.5Iterative methods to solve the equation on the graph5.1 Classical fixed point methodDenote by H L (B) the subspace of H(B) consisting of functions linear at every edge e.Consider the mapping L: H L (B) ∩ B R → H L (B) such that Lp is a solution of the linear problem on the graph: . To get this relation, one should solve problem (4.1) for all α on the cross-section of each tube.Instead, we provided the direct COMSOL 3D simulation of the non-Newtonian flow through a cylindric tube of the length 0.1 m and diameter 0.01 m for a set of values of α.On the lateral boundary of the cylinder, the no-slip boundary condition was set, while on the inflow and outflow the normal velocity corresponding to the given flux was defined.We found out that F σ (e) ε (α) can be approximated as follows: F σ (α) = 3.3372 • 10 −10 • α 1.4164 m 3 /s.

Figure 3 .
Figure 3. Level surfaces of the velocity in a tube.

Figure 4 .
Figure 4. Curve flux/pressure slope.Continuous line corresponds to the interpolation, while the dashed line is defined by the least square method.

Figure 5 .
Figure 5. First graph example with 6 nodes/vertices O i and 6 edges e j .

Table 1 .
Values of graph points O i for the first example.The "Exact" line gives the solution found with MATLAB nonlinear equations solver Fsolve, while others show results acquired with the proposed iterative method.Here, equations were scaled by 10 10 and for the optimal convergence rate, the step size in artificial time was taken equal to 1/22.Second graph example with 23 nodes/vertices O i and 24 edges e j .

Table 2 .
Values of graph points O i for the second example.The "Exact" column gives the solution found with MATLAB nonlinear equations solver Fsolve, while others present results acquired with the proposed iterative method.Here, equations were scaled by 10 19 and the step size in artificial time was taken equal to 1/58.