A mathematical model of an arterial bifurcation Текст научной статьи по специальности «Математика»

URAL MATHEMATICAL JOURNAL, Vol. 5, No. 1, 2019, pp. 109-126

DOI: 10.15826/umj.2019.1.011

A MATHEMATICAL MODEL OF AN ARTERIAL BIFURCATION

German L. Zavorokhin

St. Petersburg Department of the Steklov Mathematical Institute, Russian Academy of Sciences 27, Fontanka, St.Petersburg, 191023, Russia [email protected]

Abstract: An asymptotic model of an arterial bifurcation is presented. We propose a simple approximate method of calculation of the pressure drop matrix. The entries of this matrix are included in the modified transmission conditions, which were introduced earlier by Kozlov and Nazarov, and which give better approximation of 3D flow by 1D flow near a bifurcation of an artery as compared to the classical Kirchhoff conditions. The present modeling takes into account the heuristic Murrey cubic law.

Keywords: Stokes' flow, Bifurcation of a blood vessel, Modified Kirchhoff conditions, Pressure drop matrix, Murrey's law.

Introduction

The circulatory blood system in the human body contains about a billion bifurcation nodes, which help to distribute blood throughout the body and return to the heart. Most of the nodes belong to the capillary system, the smallest blood vessels, which are peculiar to their functioning laws, cf. [11]. Arterial and venous systems have also a considerable large number of nodes. The arterial nodes, in which the blood flows under significant pressure, are the most vulnerable elements of the circulatory system. Changes in the vessel wall geometry cause stress concentrations, provoking local damage, which in turn could create an obstacle to blood flow.

Our goal is to propose a simple approximate method of calculation the pressure drop matrix Q [18]. This is achieved by using classical shape optimization techniques. The pressure drop matrix was introduced in [12, 18] as an integral characteristic of a junction of several pipes with absolutely rigid walls. It appears that the elements of this matrix are included in the modified Kirchhoff transmission conditions, which describe more adequately the total pressure loss at the bifurcation point of the flow passed through the corresponding junction of the pipes, see [3, 9, 10].

In [10], a one-dimensional model of Stokes' flow at a junction of thin vessels with rigid walls for a fixed flow of the fluid at the inlet cross-sections and fixed peripheral pressure at the outlet cross-sections was developed. With the help of the pressure drop matrix Q apart the Neumann conditions (given flux) and the Dirichlet conditions (given pressure) at the external vertices, the one-dimensional Reynolds equations on the edges of the graph are supplied with transmission conditions at the internal vertices containing a small parameter h, where h is the thickness of the vessels, and passing into the classical Kirchhoff transmission conditions as h ^ +0. It has been established that the pre-limit transmission conditions ensure an exponentially small relative error O(e-p/h), p > 0, in the calculation of the solution to the three-dimensional problem, but the limit Kirchhoff conditions only give polynomially small relative error O(h) for the pressure and O(h3) for the velocity vector.

For obvious reasons (cf. the discussion in [9]), the bifurcation node functioning is predetermined by two of its characteristics, that is, the lumen shape and the elastic wall properties. Moreover,

the influence of each of them causes distortion for the classical Kirchhoff conditions (that are only valid in the first approximation) in a one-dimensional model. The Kirchhoff transmission conditions require the vanishing of the total flux and the pressure matching at the adjacent ends of the vessels. So shrinking/dilatation of walls violates the first Kirchhoff law as a result of provoking the local currents. This aspect of the node functioning remains outside the scope of this paper. At the same time, it is easy to take into account the uniform shrinking/dilatation of the lumen of the elastic node (the change of coordinates x ^ hx transforms the pressure drop matrix Q ^ h—3Q, see the discussion in [9]) and therefore, neglecting elastic side effects, we restrict the study to a connected channel system with rigid walls. Of course, the problem considered finds applications in other matters of fluid mechanics, not only in hemodynamics.

In Section 1, we consider the Stokes system in an unbounded domain with cylindrical outlets to infinity (see, e.g., [1, 2, 4, 14, 21]) and prove the unique solvability of the problem stated in the class of "energy" solutions, see [12] for details. For obtaining the asymptotic behavior of the solution, we exploit special homogeneous solutions to the Stokes problem with non-zero flux and with a linear growth in the pressure at infinity (see [18]). As a consequence, we obtain a definition of the symmetric pressure drop matrix Q, which plays a crucial role in the functioning of the bifurcation node.

In Section 2, we apply asymptotic analysis for the study of the bifurcation of thin vessels with rigid walls. Using the method of matched asymptotic expansions (see, e.g., [6, 22]), the modified Kirchhoff transmission conditions have been derived.

Section 3 refers to some basics for implementation in computer modeling. The introduced matrix Q of pressure drops depends on the shape of a bifurcation region and possibly on its elastic properties, see Assertion 1 in Section 1. To find such dependence for real vessels it is possible only in solving inverse problems based on data obtained by the MRI method. Nevertheless, we propose an elementary procedure of finding these coefficients in the case of a three-dimensional Stokes system (cf. [9] in the case of a two-dimensional Stokes system).

The heuristic Murray law is used to describe the radii of bifurcating vessels (see [16] and [5, Ch. 3.3]). By minimizing the cost function, which is the sum of the friction power loss and the metabolic energy proportional to blood volume, Murray derived an optimal condition for a vascular bifurcation. Murray's law states that the cube of the radius of a parent vessel equals the sum of the cubes of the radii of the daughters, namely,

R0 = R+ + R—. (0.1)

It means that the radius of at least one of the two vessels R± after bifurcation differs a little from the radius R0 of the vessel before the bifurcation. This law was verified experimentally [7, 8] but it did not get any theoretical justification. However, a number of optimization problems in branching arteries were solved by using this law (see, e.g., [5, Ch. 3.3]). We emphasize that the third powers of the radii often appears in haemodynamics.

Murray's law was a reason for several attempts (see [5, Ch. 3.3]) to derive a dynamic law of distribution of the blood flux in a bifurcating artery but should not be considered as successful. Such local law cannot exist since the flow leaving the main vessel is always defined by global data. The most frequently used in applied models of arterial flows are the classical Kirchhoff transmission conditions and their interpretations (see [17] on this topic and reference therein, and some reviews in [5, 20]). In any case, such conditions take into account neither the elastic properties of a vessel nor the most of the geometric characteristics of the bifurcation region. The modified Kirchhoff conditions have no these lacks.

In Section 4 it will be explained why the modification of the second Kirchhoff law by means of the pressure drop matrix unexpectedly deeply increases the accuracy of approach for three-dimensional fluid flow in a system of thin channels by the one-dimensional Reynolds-Poiseuille

model. In addition, we present a simple relation for a distributional law of Murrey's type.

1. Statement of the problem

1.1. Domains with cylindrical outlets to infinity and functional spaces

We introduce the domain Q with three cylindrical outlets to infinity (see Fig. 1). Let Q be an open unbounded domain with Lipschitz boundary dQ admitting the representation

Q = Q' U Q0 U Q+ U Q-, where Qa n Q^ = 0 for a = = 0, ±.

Here, Qa = {xa = (ya, za) : ya G ua, za > La} in a certain Cartesian coordinate system xa = (ya ,za) in R3, where ya are the variables in the cross-section of the outlet Qa, za is the variable along the axis of Qa, and is a bounded domain in R2. The bounded domain Q' is given by Q' = {x G Q : za < L} for certain L, L > maxa La. Henceforth, x = (xi, x2,x3) is a global coordinate system in R3 related to the whole domain Q.

x,

Figure 1. Artery bifurcation (the domain Q). We define L2,^ (Q) as the space of measurable functions in Q with a finite norm

u||L ß(O) =(f |u(x)|2dx + £ Í |za|2ß|u(ya,za)|2dyadza)1/2. VJo! ."Ti Joa J

a=0,±

If ^ = 0, we will use the usual notation L2 (Q) for this space.

By using the Sobolev space H1 (Q) together with L2j1 (Q), we introduce the space of real-valued vector-functions in Q

H(fi) = {u = (ui,U2,U3) G (H1 (fi))31 divu G L2,1 (fi)}

with the norm

u

I H(O)

(|Vu(x)|2 + |u(x)|2 )dx + V / |za |2 |div u(ya ,za )|2 dya dza.

a=0,±

Let also H0 (Q) be the subspace in H(Q) consisting of vector-functions equal to zero on dQ. The dual space of H0 (Q) is denoted by (H0 (Q))*.

z

2

O

1.2. Formulation of the problem

Consider the Dirichlet problem for the stationary Stokes system with nonzero divergence

—vAu(x) + Vp(x) = F(x), -divu(x) = G(x), x € Q, (1.1)

u(x) = 0, x € dQ. (1.2)

Here, u(x) = (ui(x), u2(x), u3(x)) is the velocity field, p(x) is the pressure, and v > 0 is the viscosity of the fluid, which is assumed to be constant.

In order to define a weak solution of problem (1.1), (1.2), we introduce a bilinear form on H(Q):

3 f

i(u, w) = 2_] / VujVwjdx. j=i ^

So, if (u,p) is a classical solution of (1.1), (1.2), then multiplying the first equation in (1.1) by w € H0(Q) and integrating over Q, we obtain

va(u,w) — p div w dx = Fw dx for any w €H0(Q). (1.3)

Jo Jo

A weak solution of problem (1.1), (1.2) is a pair (u,p) € H0(Q) x L2;-1(Q) satisfying the integral identity (1.3) for all w € H0(Q) and the equation —div u = G in Q, where F € (H0(Q))* and G € L2;1(Q) are given.

To prove the main result of this section, we need the following statement.

Lemma 1. For arbitrary g € L2;1(Q) subject to

/ g(x)dx = 0, (1.4)

jo

there exists a vector-function u € H0(Q) such that —div u = g in Q and

llullH(o) < c\\g\\h2,i(o). Here, c is a constant independent of g.

Proof. We first consider an auxiliary problem; namely, we look for a solution of the boundary value problem

—div ua(xa) = na(ya)Ga (za), xa € Qa, (1.5)

ua(xa) = 0, xa € dQa, (1.6)

where

na(ya) € C0°>a), f na(ya)dya = 1 and Ga(za) =f g(ya,za)dya.

J Wa J Wa

One can verify directly that the vector-function

f <x

ua(xa) = (0, 0, na (ya)wa(za)), wa(za) = — Ga(t)dt,

Jza

solves problem (1.5), (1.6). Moreover, by Hardy's inequality (see [13, Theorem 5.2])

l|wa||L2(La>~) < W \za|2|Ga|2dza < C||g||L2,1(La>~). (1.7)

J La

Using (1.7), we obtain

||ua\|H(0a) < ca\\g||L2,i(O«). We are looking for a solution to —div u = g in the form

u = £ Xa(za )ua(ya ,za)+ U (x) (1.8)

a=0,±

for the equation —divu = g in Q. Here, the cut-off functions xa are defined by Xa(za) = 1 for za > 2La, Xa(za) = 0 for za < La. The function U in (1.8) satisfies the equation

—divU = g + div( Xa(za)ua(ya,za)^ = g € L2,i(Q)

a=0,±

in Q. It is easy to see that

/ gdya = 0 for za >La. (1.9)

J Wa

Let us introduce a local covering of Q. Let

Qa = {(ya, za) : ya € Ua, La + j — 1 < za < La + j + 3/2}, j = 1, . . .

Then

TO

Q = Q' + ££Qi.

a j=1

Let us consider the partition of unity corresponding to this covering:

to

1 = ('(x) + £ £ ^(za), where j € C^(La + j — 1, La + j + 3/2) a j=i

and <// is a smooth function supported in Q1. We can write

g = g' + 52 ga, where gja = (jg, g' = ('(x)g. a j=i

By (1.9) and (1.4), we have

/ gfdxa = 0, and, / g'dx = 0. Jo« jo'

So, we obtain the problem in each bounded domain Qa:

—div Ua = ga in Qa, (1.10)

Ua = 0 on dQa, (1.11)

and a similar problem for U' in Q'. These problems have solutions Uj € H0(Qa) (see [21]) which satisfy

||Uja||H(oa) < C \ |gj* \ \l2,i(o«) .

Thus, the vector-function

ro

U = U' + j €Hq(H)

a j= 1

solves problem (1.10), (1.11) and satisfies

H(Q) < Cl( ||U '||H(q' ) + EE llj Ikj < C2||g||L2il (Q).

a j=1

The proof is complete. □

Now we prove the existence and uniqueness of weak solutions to the boundary value problem (1.1)-(1.2), cf. [12].

Theorem 1. Suppose that F e (%Q(H))* and G e L2;1(H) is such that

I G(x)dx = 0. (1.12)

Jq

Then there exists a uniquely (up to an additive constant in the pressure p) determined weak solution (u,p) e hq(Q) x L2;-1(H) of problem (1.1), (1.2). Furthermore,

||u||H(Q) + ||P||L2,-1 (Q) < c (||F 11 (%o(Q))* + ||G||L2,1 (Q)) . Here, c is a constant independent of F and G.

Proof. Let V(x) = u(x) — w(x), x e H, where u(x) is a solution to the problem (1.1), (1.2) and the conditions of Lemma 1 for w(x) are satisfied, namely,

—div w(x) = g(x), x e H, w(x) =0, x e dH.

Then, we arrive at the Stokes problem

—vAV(x) + Vp(x) = F —div V(x) =0, x e H, (1.13)

V(x) = 0, x e d H, (1.14)

where (HQ(H))* 9 F(x) = F(x) + vAw(x). The weak formulation of problem (1.13), (1.14) is as follows: find V e (H) such that

va(V, W)= f F dx for any W e HQiv(H). (1.15)

Q

Here, Hdiv (H) = {W e HQ (H) : div W = 0 in H} . By the Riesz theorem, there exists a uniquely determined vector-function V eHdiv (H) such that ( ) is satisfied and

IIV||H(q) < c||fP||(Ho(q))*.

Here, we introduce

l2,i(H) = {g(x) e l2,i(H) : Jq g(x)dx = 0}.

ro

By Lemma 1, for any q € L2)1 (Q), there exists a vector-function uq € H0(Q) such that —div uq = q in Q and

\|uq\\h(o) < c\\q||L2,i(o).

We consider the functional

G(q) = / Fuq dx — va(V,uq) (1.16)

o

on L2,1(Q). In virtue of

|G(q)| < c( ||F ||(Ho (o))* + ||V ||„(o)) l l uq | | H(o) < c | | F | | (Ho (o))* l l q 11 L2,1 (o),

the linear functional G(q) is continuous on L2;1(Q) and, by the Riesz theorem, there exist a unique element p of the dual space L2;-1(Q) such that

G(q) = pq dx for all q € L2;1(Q).

o

Then, by (1.16), we get

va(V, uq) + pq dx = / Fuq dx,

oo

or

va(V,uq) — p divuq dx = Fuq dx for all uq € H0(Q).

oo

We have the estimate for

iipiil2-i (n) = suP

L2,i(H): ||q||L21(Q) = 1

/ pq dx n

sup

L2 ,i(n): ||q||l2 1(n) = 1

/ Fuq dx — va(V, uq) n

< c||F||(Ho(n))*

Finally, then (u,p) = (V + w,p) is the required weak solution, and the estimate

\|u| \h(o) + \\p||L2 ,-i(o) < c (||F \ \ (Ho(o))* + ||G||L2 , 1 (o)) is fulfilled. □

Remark 1. Consider a non-homogeneous Dirichlet problem for Stokes system, i.e., equations (1.1) are supplied with the boundary condition

u(x) = H, x € dQ, (1.17)

where H € H(Q), and, instead of (1.4), we require

/ G(x)dx + H(x) ■ nd r = 0, Jo Jdo

where n is the unit outward normal to dQ. Substituting u(x) = v(x) + H(x) into (1.1), (1.17), we obtain

—vAv(x) + Vp(x) = f (x), —div v(x) = g(x), x € Q, v(x) = 0, x € d Q,

where f (x) = F(x) + vAH(x) € (%0(Q))* and g(x) = G(x) + div H(x) € L^Q) verifies (1.4). Now, the application of the previous theorem gives the existence of a pair (v,p) € H0(Q) x L2;-1(Q) solving problem (1.1), (1.17) and satisfying the estimate

\|v| \h(o) + \\p||L2 ,-i(o) < ^||f \ \(Ho(o))* + \\g||L2 , 1 (o) + ||H\\h(o^ . Moreover, p is defined up to an additive constant.

1.3. Asymptotics of the variational solution

Let the right-hand sides in (1.1), (1.2) satisfy

I |F(x)|2dx + V/ |F(xa)|2e2az"dxa < to (1.18) JQ' ^ JQa

and

f |G(x)|2dx + V/ |G(xa)|2e2az"dxa < to, (1.19)

JQ' a JQa

where a is a positive number. Let also G be subject to (1.12). Then, according to Theorem 1, problem (1.1), (1.2) has a solution (u,p) € H0(Q) x L2;-1(Q). We can conclude that this solution satisfies the following asymptotic representation at infinity:

(u,p)= Y, XaCa (0, 1)+(Vp), (1.20)

a=0,±

where xa = Xa(za) are smooth functions equal to 1 for za > La + 1 and 0 for za < La, (V,p) are exponentially decaying terms as za ^ to, and ca are real constants. Since this solution is defined up to an additive constant, we can (and will) assume c0 = 0. Then the solution is unique, see Remark 1.

The remaining part of this section is devoted to finding formulas for evaluation of the constants c+ and c_. For this purpose, we need solutions of homogeneous problem (1.1), (1.2), which have a linear growth at infinity. More precisely, we introduce two linear independent solutions (V±,P±) which have the following asymptotic representations (see [18]):

(V± ,P±) = —Xo(V0, P0) + x±(V±, P ±) + (v± ,p±), (1.21)

where (Va, Pa) is the Poiseuille flow in the cylinder Qa, i.e., VO* (x) =0, i = 1,2, Pa(x) = —A-1za, and Vaa = Vaa(ya) solves the following Dirichlet problem in wa:

AVO* = —A-1 in Wa, VZa = 0 on dwa. (1.22)

The normalizing constant ca is choosing to satisfy

/ vza(ya)dya = 1. (1.23)

Jua

In the most important case of the circular cylinder, i.e., wa = {ya : |ya| < ra}, we have

2(r2 — |ya I2) —8V

V£(*) = J > = (1.24)

" ' a " ' a

The remainder term (v± ,p±) in (1.21) satisfies the problem

—vAv± + Vp± = f ±, —divv± = g± in Q, v± = 0 on d Q,

where the right-hand sides

f± := vA(x0V0 — x±V±) — V(X0P0 — X±P±) (1.25)

and

g± :=div(x0V0 — X±V±) (1.26)

have compact supports. To verify condition (1.2) in Theorem 1 for g±, we apply the Gauss theorem for the domain Qr = {x € Q : za < R}, where R is a sufficiently large number, and obtain

f g±dx = lim f div(xoV0 - x±V±)dx =/ Vz°o(y°)d£ -/ V±± (y±)dE = 0.

Jn R^J Qr o J„± z

Here we used the normalization condition (1.23). Therefore, (v±,p±) admits the asymptotic representation (1.20), where c° = 0 and (V,p) exponentially tends to zero as za — to. Now we can present formulas for calculation of coefficients in (1.20)

Theorem 2. Let the functions F and G satisfy (1.18)-(1.19) and let the asymptotic formula (1.20) be valid with c° = 0. Then

c± = / (FV± + GP±) dx. (1.27)

Jn

Proof. Let Qr be the same domain as before. Multiplying equations (1.1), (1.2) by (V±, P±), integrating over QR and using Green's formula, we obtain

/ ((-vAv + Vp) V± - divvP±) dx

JnR

= £ / ( - v (V±dza vza - ^ Vz± ) + (pVz± - P± V^a ) )

a •'"a

Taking here limit and using asymptotic formulas for (V±,P±) and (1.20) we arrive at (1.27). □

z*=R y

Applying formula (1.27) to the solution (v±,p±) of problem (1.1), (1.2) with the right hand sides given by (1.25), (1.26), we obtain the representations

(v± ,p± ) = X±Q±±(0,1) + xtq±t(0, 1) + (v± ,P± ), (1.28)

where the coefficients are evaluated according to the formula

qyt = / (fYVT + g1 PT) dx, y,t = ±. (1.29)

Jn

From (1.21) and (1.28), we get the following representations:

(V ±, P ±) = -xo(v0, P0) + x±(v ±, P ± )

+X±Q±±(0,1) + xtq±t(0, 1) + (v±,P±) (1.30)

with the remainders (v±,p±) exponentially decaying at infinity. Note that a straightforward calculation gives the equality QYT = QTY. The coefficients QYT in expansion (1.30) of the pressure at infinity in Q± form the symmetric (2 x 2)-matrix Q called the pressure drop matrix. Another approach to introducing the matrix Q was presented in [18].

Assertion 1. Since the replacement p — p = v-lp eliminates the viscosity v in the Stokes system (1.1), the matrix Q admits the representation Q = vQ, where Q depends only on the geometry of the set Q (the bifurcation node).

Remark 2. Solution (1.30) describes the flow of the fluid (blood) with unit flow through the cross-section w0 in the direction of the bifurcation node and its outflow through the cross-section w±, that is, the flow through the other cross-section is zero. The constant component in the representation of the pressure at infinity in the supplying blood vessel Q0 is also made to be zero, but P± contains the linear component A-^z0, which together with — A^.1 z± ensures a pressure drop between the infinitely distant cross-sections in the blood vessels Q0 and Q±. The summands with the coefficients Q7T in (1.30) also make their contribution, however, it is finite and, therefore, the (2 x 2)-matrix Q = (Q7T)Y,T=± composed of them is called the pressure drop matrix. We note that the absence of the summand x0Q0± in (1.30) fixes the pressure P± and, therefore, the matrix Q itself, while earlier the arbitrariness in the solution to problem (1.1), (1.2) was eliminated by the condition c0 = 0. Due to conditions (1.18), (1.19), we apply the fact that some external actions F and G in the model are negligible on cylindrical outputs at infinity, since we study in principal the problem with the bifurcation of an artery.

2. Asymptotic analysis of the bifurcation of thin channels with rigid walls

Let the bifurcation be characterized by a common small parameter h > 0. More precisely, suppose that

Qh = Q' U Qh U Q+ U Q-,

Q£ = {x : |ya| <hr„, <U , a = 0, ±,

ra > 0 and > 0 are certain fixed radii and lengths, respectively, the radii r0 and r± are comparable, and h << max{1, r-V±}.

In the domain Qh, consider a stationary Stokes flow. The velocity vector uh and the pressure ph satisfy the Stokes system

-vAxuh(x) + VxPh(x) = 0, -div uh(x) =0, x € Qh, (2.1)

with the no-slip conditions

uh(x) = 0, x € Xh = d Qh n d Q. (2.2)

on the lateral boundary.

Nonhomogeneous boundary conditions at the ends = {x € d Qh : = 1a}, describing the inflow and outflow of the fluid, will not be used, but we can take them the same as in Section 4. Asymptotic ansatzes for a flow in a pipe Q^ (in what follows, it will be called the Reynolds-Poiseuille ansatz) have the form

ph (x)= p«(za ) + ..., (2.3)

h2

(x) = 0 + • • • , i = 1,2, Uhza(x) = — {rl - h-2\ya\2)dzapa(za) + • • • (2.4)

Here, and uh« are the velocities along the axes yf and The flux through the cross-section in the positive direction is evaluated as

I lit (X)\2~=2dya = h4^ =: h4Aa. (2.5)

|ya|<hra

The unknown functions a = 0, ±, satisfy the classical one-dimensional Reynolds equations

-Aadz2.p«(za) =0, € (0, la). (2.6)

Boundary conditions at the endpoints za = la are not needed, but to adjust the solutions pa(za) at the common point O with coordinates za = 0, we apply the method of matched asymptotic expansions (see, e.g., [6, 22]).

Theorem 3. The solutions pa, a = 0, ±, of the Reynolds Eqs. (2.6) satisfy the modified Kirchhoff transmission conditions at the bifurcation point za = 0

po(0) = p±(0) + h Y, Q±tAt0ztpT(0), (2.7)

T = ±

£ Aadza Pa (0) = 0. (2.8)

a=0,±

Proof. We stretch the coordinates with respect to the center of the bifurcation node:

x £ = h-1(x -O). (2.9)

Taking the formal limit at h = 0, we transform Qh into the unbounded domain Q1 with three outlets to infinity having the shape of the semicylinders Q1a = {£ : |na| < ra,(a > 0}, where a = 0, ±, and £a = (na,Za) are local coordinates obtained by stretching the coordinates (ya,za) (cf. (2.9)).

As usual, for the application of the method of matched asymptotic expansions, it is required to find all solutions to the Dirichlet problem for the Stokes system in Q1, which may grow at infinity no faster than linear functions. One of such solutions is obvious, that is, the constant pressure

P0(£ ) = 1, V0(£)=0.

There are two more solutions (1.30) generated by the unit flux of fluid in Q10 from infinity, which is compensated by an equal flux in Q1± at infinity as well. Now we will use the procedure of matching inner and outer expansions. As an inner expansion, we take the linear combinations

h£>±V ±(£), (2.10)

±

ah0P0(0 + E a±P±(£). (2.11)

±

The additional factor h is used in the velocity expansion (2.10) according to the general algorithms of construction of boundary layers (see, e.g., [15, Ch. 4]); it is needed, because the vector uh in (2.1) is differentiated twice, but the scalar ph only once.

By using the Taylor formula with respect to za, we rewrite the outer expansions (2.3) and (2.4) in the form

ph(x) = pa (0) + zadzapa (0) + ... = Pa (0) + h(a0zapa(0) + . . . , (2.12)

h2 \v I2 h2

vh{x) = _(r2 _ nM)dzapa{Q) + ... = _(r2 _ [r]«\i)dzapa{0) + ... (2.13)

Comparing the coefficients of constant and linear functions in representations (2.12) and (2.11), we take into account (1.21)-(1.24) and (1.30) and arrive at the relations

po(0)= ah, p± (0) = ah + Q±+a+ + Q±-a-,

hAodza po(0) = a+ + a+, hA± dz* p± (0) = -a±,

where the quantities Aa from (2.5) and (2.6) and (1.21)-(1.24) are used. These relations guarantee the matching of the leading terms in the asymptotics (2.4), (2.10), and (2.13) and imply the relations

ah = Po(0), a± = -hA±dza p±(0)

and (2.7), (2.8). □

If we logically take the limit as h — 0 in relations (2.7) and (2.8), then the latter terms in (2.7) vanish and we obtain the classical Kirchhoff transmission conditions

P±(0) = P0(0), E Aadza pa(0) = 0,

a=0,±

which mean the continuity of the pressure and the vanishing of the total flux from the bifurcation point (cf. (2.5) and (1.23)).

Let us explain why it is preferable to keep the last two terms in (2.7) despite they are small, see the discussions in [9, 10]. The solutions to the differential equations (2.6) are linear functions

pa(za) = cP(h) + zaci(h), a = 0, ±. (2.14)

The coefficients c^, n = 1,2, can be found from the Kirchhoff conditions (2.7), (2.8), and boundary conditions at za = see Section 4. The quantities c^ = c^(h), n = 1,2, are rational functions of the small parameter h, since they are determined from a system of linear algebraic equations with invertible matrix, polynomially (linear) dependent on h. The functions

ph(x) = pa(za) = c£(h) + zaci(h),

h2

u:;?{x) = 0, i = 1,2, uhza(x) = — (r2 - h-'2\ya\2)cla{h)

constructed according to (2.12), (2.13), and (2.14) form the Poiseuille flow in the cylindric parts of the vessels Q^ and, hence, they exactly satisfy the Stokes system (2.1) and the no-slip condition (2.2) restricted to the corresponding parts of the literal surface.

Near the bifurcation region of Qh, there appears a boundary layer (we go over from the method of matched asymptotic expansions to the method of compound asymptotic expansions; cf. [15, Ch. 2]). By definition, relations (2.7) and (2.8) become conditions of an exponential decay for the boundary layer. As a result, we derive that the constructed approximate solution leaves in the problem (2.1), (2.2) exponentially small discrepancies as h — 0. Thus, using the more complicated coupling conditions (2.7), (2.8), we obtain estimates of the asymptotic remainders with majorants Ce-p/h, p > 0, whereas the application of the Kirchhoff conditions delivers remainders of order O(h) in the asymptotic representation of pressure and of order O(h3) for the velocity.

3. Asymptotics of the pressure drop matrix

3.1. Elementary procedure for finding the pressure drop matrix. Murrey's law

Consider an infinite three-dimensional symmetric channel Y with a bifurcation, which is depicted in Fig. 2, where its geometric parameters are shown. The walls are assumed to be rigid and the angle 0 is small. The vertical dot-and-dash divide Y into four parts: three semi-infinite cylinders Q± and Q0 and a middle part Q^, which is located between the points z = 0 and z = — L in the figure. The length and width of the middle section Q^ are evaluated as follows:

L = (sin60-1(2R — R0cos0) « (sin60-1(2R — R0), R = R±, ( )

2a^ (z)=2(R0 — z sin 0).

Figure 2. Symmetric bifurcation of a three-dimensional channel.

We assume that the cross-section of the middle part ^ is an ellipse whose semi-axes are connected by the relation

b.(z) = (1 + Pz)a.{z), [—L,0], with f3 = L~1{2 + V2)-1.

Let us find approximate formulas for solutions (1.30) of problem (1.1), (1.2) in the three-dimensional junction Y. Its symmetry with respect to the horizontal axis allows us to consider only one of the solutions, e.g., (V,P) = (V + ,P +).

Theorem 4. Let the radii Ra, a = 0, ±, of the supplying and accepting channels obey Murrey's law (0.1). Then the pressure drop matrix Q in (1.30) is positive definite.

Proof. Owing to the large parameter (sin6)-1, we can use a three-dimensional Reynolds-Poiseuille ansatz (2.3), (2.4) ignoring the boundary layer effect, which brings an error of order O(1), small with respect to O((sin 6)-1). Thus, we seek the one-dimensional distribution of the pressure in the form

P(x) « P°(z) = A-1 z for z e (0,

P(x) « P + (z) = A-1 z + Q++ for z e (-to, -L.), (3 2)

P(x) « P-(z) = Q+_ for z e (-to, -L.), ( . )

P(x) « P• (z) for z e (—L., 0).

In the three-dimensional case, the coefficients, cf. (2.5) and (1.21)-(1.24), are found by the formulas

tt a!(z)(l + /3z)3 ttR% _ ttE^

A'[Z)~ Av {! + {! +pzYY 8!/' 8v '

(see, e.g., [19] for calculation of an explicit formula in case of ellipse of the coefficient A., which is proportional to the torsion rigidity of a section), and the Reynolds equation on (—L., 0) takes the form

-8Z (A. (z )dz P. (z)) = 0. (3.3)

— /

Taking into account (3.2), solving the Reynolds Eq. (3.3) for unknown function P• and satisfying the Kirchhoff transmission conditions at z = 0 and z = —L., we find the quantities QYT, j,r = ±, in (1.30):

Q+_ » P'{-L.) »--^-^-(iî.iîo),

n sin ^ 4v (3.4)

Q++ » P\-L.) + A~lL. «--—-<f>++(R, Rq).

n sin U

Here, the functions have explicit but cumbersome expressions and do not depend on the small parameter sin 9.

Using Sylvester's criterion for the pressure drop matrix Q, we state the assertion of Theorem 4.D

z=0 z=-l.

Figure 3. Asymmetric bifurcation of a three-dimensional channel.

A three-dimensional channel Y with branching is depicted in Fig. 3, where all geometric parameters can be found. We accept the notation and assumptions from the symmetric bifurcation. Formulas (3.1) are replaced by

L. « (sin 9)-12R, 2R(z) = 2Ro - z sin 9. (3.5)

Here we assume that the cross-section of the middle part ^ is a disc of radius R^ (z), z G [—L, 0]. Since there is no symmetry now, we must construct both solutions (V±, P±).

Proposition 1. The same Theorem 4 is valid for the asymmetric bifurcation (see Fig. 3).

Proof. We start with (V, P) = (V +, P +), the flux of which goes to the left channel of width 2R0 and, similarly to (3.2), we obtain

P(x) « P0(z) = A-1 z for z G (0, +to),

10

P(x) « P+ (z) = A-1 z + Q++ for z G (-TO, -L.),

(3.6)

P(x) « P-(z) = Q+- for z G (-TO, -L.),

P(x) « P• (z) for z G (-L., 0).

We have the relations

Equation (3.3) is valid, and the Kirchhoff transmission conditions at z = 0 and z = -L. show that

Sv f 1___1_\ 8u f 1___1_ 6R\

g+-~37rsind\(Ro + R)3 R3)' Q++~ 3nSme\(R0 + R)3 R30 + R40)' [ '

For the solution (V, P) = (V-, P- ), the flux of which goes to the lower (see Fig. 3) channel of width 2R, the two middle relations in (3.6) must be changed according to

P (x) « P- (z ) = A-1 z + for z G (-TO, -L. ),

P (x) « P+ (z ) = Q-+ for z G (-to, -L. ).

In this case, we obtain

q_I__n _I__J- + -ËA f38)

3vrsin6\(R0 + R)3 R3)' 37rsin0 \(R0 + R)3 R30 R3 J ' '

Using Sylvester's criterion, we state that the pressure drop matrix Q is positive definite. □

Remark 1. A possibility to obtain approximate formulas (3.4) and (3.7), (3.8) for the pressure drop matrix Q is provided by the fact that, for small 6, the junction elongates and has length O((sin6)-1) (see expressions (3.1) and (3.5) for L.). At the same time, the limit passage 6 — +0 fails, i.e., one can use these results only for "not very small" angles 6.

Interesting to note that the entries (3.4) and (3.7), (3.8) of the matrix of pressure drops are inversely proportional to the cube of the radii in accordance with Murrey's distributional law (0.1).

4. Pressure drop matrix and modified Kirchhoff transmission conditions

Let us truncate cylindrical outlets in Q and assume

Qh = Q' U Qh U Q+ U Q-, Qa = {x : lyal < ra, za <h~1la] , a = 0, ±, (4.1)

where h > 0 is a small dimensionless parameter, and ra > 0 and la > 0 are certain fixed radii and lengths, respectively. In the domain Qh, we define the homogeneous (F = 0, G = 0) Stokes equations (1.1) and, on its lateral surface £h = dQh n dQ, we impose the homogeneous (H = 0) no-slip conditions (1.2) (hereinafter, we refer to these relations implying that they are restricted to these sets). On the truncated surfaces ra = {x : |ya| < ra, za = h-1 la} , assign the following conditions:

4(x) = 0, i = 1,2, uhz0(x) = -*h(y°), x e rh, (4.2)

vt (x) = 0, i = 1,2, -vdzTuhzT (x)+ ph(x) = p~, x e rh, t = ±. (4.3)

Here, is the Prandl function, that is, the solution of the Dirichlet problem for the Poisson equation (1.22). In other words, at the inlet cross-section of the vessel Qh, the incoming unit flux of fluid is assigned and, on the allocated ends of the outlet cross-sections of the vessels Q±, peripheral pressure p^ is set. At the same time, the compression of coordinates by h-1 times transforms the problem stated to the usual problem of the blood flow through the bifurcation node of thin vessels, which walls, as already explained, it is assumed to be rigid (see Section 2, cf. [9]). In the new coordinates, the vessels become the smaller radii hra and the fixed lengths la. We emphasize that the problem (1.1), (1.2), (4.1), (4.2), (4.3) is still included in the symmetric Green formula in Qh. Its interpretation in the framework of the weighted spaces technique with detached asymptotics is given in [18]. The conventional Reynolds-Poiseuille ansatzes (see Section 2)

uh (x) = —(8v )-1 n*a(ya )eza dza Pa(za) + ..., ph (x) = Pa(za) + ... (4.4)

(where eza is the unit vector of the axis za directed away from the node) on the vessels Qa, after the substitution into (1.1), (1.2), (4.2), (4.3), generate the differential equations and boundary conditions

2«Pa = 0, for € (0,h-11a), A« = ^v)-1^,

0 _ . _ ^^ _ u-1i

—A°3zop° = 1 at z° = h 1l°, p± = p^ at z± = h 1l±

for unknown functions pa, a = 0, ±. The latter functions linearly depend on the longitudinal variable za, and, thus, it is possible to remove the dots replacing lower-order asymptotic terms in the right-hand sides of Eqs. (4.4).

As an approximate solution of the problem stated in Qh, we take the sums

uh = ha+V + + ha-Vph = a+P + + a-P- + ah, (4.5)

where (V±,P±) are introduced special solutions (1.21) and the latter term refers to the constant pressure. Using the method of matched asymptotic expansions in the interpretation of ansatzes (4.4) as the outer expansions and the linear combination (4.5) as the inner expansions (see Section 2, cf. [9]), we satisfy the boundary conditions (4.2), (4.3) up to exponentially small terms as h — +0 and get the following relations:

1= a+ + a-, (4.6)

= ah - h_1Lrah + ^ Qraaha, t = ±, (4.7)

a=±

8v

where Ca = —jla, a = 0, ±, and henceforth C = diag{£+, £-}.

Let e = (1,1) and ah = (a+, a-) be columns. In virtue of (4.7), we deduce

- ah)e = (Q - h-1L)ah,

hence,

ah = (h-1L — Q)-1e(ah - p~), and thus equality (4.6) rewritten in the form e ■ ah = 1 leads to the relations

1= Th(ah - ),

Th = e ■ (h-1L - Q)-1e = h-1 e ■ L-1(/ - hQL-1)-1e = h-1e ■ (L-1 + hL-1 QL-1 + O(h2))e. We finally find that

ah = + T-1,

Th = h-1(to + ht1 + O(h2)), to = e ■ L-1e > 0, ¿1 = e ■ L-1QL-1e, (4.8)

so, we get

ah = + ht-1 (1 - ht-1t1 + O(h2)). (4.9)

Thus, the pressure at the "input" rh up to the smaller terms is equal to

h-1Lo + + ht-1 - h2t-2t1 + O(h3). (4.10)

The first term of (4.10) is the pressure drop, which provides the unit flux delivery to the artery bifurcation, the second term is also positive, it is necessary to supply the fluxes to the points z± = h-11±, and the third term, the sign of which depends on the pressure drop matrix Q, corresponds to just the shape of the node.

Proposition 2. Let a± = 1/2 ±bh/2 in accordance with (4.6), and let bh € [-1,1] be the factor of flux distribution. Then the partial Murrey cubic law

= -V1"V1 (4.H)

is satisfied provided that l+/r+ = /r_.

Proof. From (4.7), we deduce

0 = -h-1L+(1 + ) + h-1L_(1 - + Q++(1 + +Q+_(1 - - Q—(1 - - Q_+(1 +

or

-h-1bh(L+ + L_) + bh(Q++ + Q— - 2Q+_) = h-1(L+ - L_) - (Q++ - Q__),

Hence,

bh= L+-L.-h(Q++-Q-)

(L+ + L_) - h(Q++ + Q— - 2Q+_)"

If we put 1+/r+ = 1_/r_, then we arrive at (4.11) as h — +0 or, similarly, the flux ratio is

1 L+-L-l + bh _ L+ + L- _ _ _ r^ 1 ~ L+-LJ ~ rt~

L+ + L_

□

5. Conclusion

For the arterial tree, under the assumption that the walls of the blood vessels are rigid, for every bifurcation node, a (2 x 2)-pressure drop matrix Q appears, and its influence on the Kirchhoff transmission conditions is taken into account. The modified Kirchhoff transmission conditions via the matrix Q depend on the geometry of the bifurcation region. We do not know the exact value of the matrix Q for any concrete bifurcation of an artery. In Section 3, we have produced the calculation scheme of an approximate determination of this matrix. With the help of Murrey's cubic law, we indicate general properties (sign definiteness) of the characteristics introduced for the certain nodes in which some information about the matrix is nevertheless accessible.

Acknowledgements

The author was supported by Linkoping University, and by RFBR grant 16-31-60112.

REFERENCES

1. Amick C.J. Steady solutions of the Navier-Stokes equations in unbounded channels and pipes. Ann. Sc. Norm. Super. Pisa CI. Sci. (4), 1977. Vol. 4, No. 3. P. 473-513. URL: http://www.numdam.org/item/?id=ASNSP_1977_4_4_3_473_0

2. Amick C.J. Properties of steady Navier-Stokes solutions for certain unbounded channels and pipes. Nonlinear Anal. Theory, Meth., Appl., 1978. Vol. 2, No. 6. P. 689-720. DOI: 10.1016/0362-546x(78)90014-7

3. Berntsson F., Karlsson M., Kozlov V. A., Nazarov S. A. A Modification to the Kirchhoff Conditions at a Bifurcation and Loss Coefficients. LiU electronic press, MAI, LiTH-MAT-R-2018/05-SE, 2018. P. 1-10. URL: https://www.diva-portal.org/smash/get/diva2:1204214/FULLTEXT01.pdf

4. Bogovskii M. E. On solution of certain problems of vector analysis associated with operators div and grad. Tr. Semin. im. S.L. Soboleva, 1980. Vol. 1, Novosibirsk. P. 5-40.

5. Fung Y. C. Biomechanics. Circulation. 2-nd ed. New York: Springer-Verlag, 2011. 572 p. DOI: 10.1007/978-1-4757-2696-1

6. Il'in A.M. Matching of Asymptotic Expansions of Solutions of Boundary Value Problems. Transl. Math. Monogr., Vol. 102. Am. Math. Soc., 1992. 281 p.

7. Kassab G.S. and Fung Y.-C.B. The pattern of coronary arteriolar bifurcation and the uniform shear hypothesis. Ann. Biomed. Eng., 1995. Vol. 23, No. 1. P. 13-20. DOI: 10.1007/BF02368296

8. Kassab G.S., Rider C.A., Tang N.J. and Fung Y. C. Morphometry of pig coronary arterial trees. Am. J. Physiol. Heart Circ. Physiol., 1993. Vol. 265, No. 1. P. H350-H365. DOI: 10.1152/ajpheart.1993.265.1.H350

9. Kozlov V. A., Nazarov S. A. Transmission conditions in a one-dimensional model of bifurcating arteries with elastic walls. J. Math. Sci. (N. Y.), 2017. Vol. 224, No. 1. P. 94-118. DOI: 10.1007/s10958-017-3398-0

10. Kozlov V. A., Nazarov S. A. A one-dimensional model of flow in a junction of thin channels, including arterial trees. Sb. Math., 2017. Vol. 208, No. 8. P. 1138-1186. DOI: 10.1070/SM8748

11. Kozlov V. A., Nazarov S. A., Zavorokhin G. L. A fractal graph model of capillary type systems. Complex Var. Elliptic Equ., 2018. Vol. 63, No. 7-8. P. 1044-1068. DOI: 10.1080/17476933.2017.1349117

12. Kozlov V. A., Nazarov S. A., Zavorokhin G. L. Pressure drop matrix of a bifurcation of an artery with defects. Eurasian J. Math. Comput. Appl., 2019. Vol. 7, No. 3. Accepted.

13. Kufner A. Weighted Sobolev Spaces. Ser. Teubner Texte zur Mathematik, vol. 31. Teubner, 1980. 152 p.

14. Ladyzhenskaya O. A., Solonnikov V. A. Determination of the solutions of boundary value problems for stationary Stokes and Navier-Stokes equations having an unbounded Dirichlet integral. J. Sov. Math., 1983. Vol. 21, No. 5. P. 728-761. DOI: 10.1007/BF01094437

15. Maz'ya V. G., Nazarov S. A., Plamenevskij B. A. Asymptotic Theory of Elliptic Boundary Value Problems in Singularly Perturbed Domains. Vol. II. Ser. Oper. Theory Adv. Appl., Vol. 112. Basel: Birkhauser, Verlag, 2000. DOI: 10.1007/978-3-0348-8432-7

16. Murray C. D. The physiological principle of minimum work. I. The vascular system and the cost of blood volume. Proc. Natl. Acad. Sci. USA, 1926. Vol. 12, No. 3. P. 207-214. DOI: 10.1073/pnas.12.3.207

17. Mynard J. P. and Smolich J. J. One-dimensional haemodynamic modeling and wave dynamics in the entire adult circulation. Ann. Biomed. Eng., 2015. Vol. 43, No. 6. P. 1443-1460. DOI: 10.1007/s10439-015-1313-8

18. Nazarov S. A., Pileckas K. Asymptotic conditions at infinity for the Stokes and Navier-Stokes problems in domains with cylindrical outlets to infinity. In: Advances in Fluid Dynamics, Quaderni di Matematica, ed. P. Maremonti, Vol. 4., 1999. P. 141-243.

19. Polya G., Szego G. Isoperimetric Inequalities in Mathematical Physics. Ser. Ann. of Math. Stud., vol. 27. Princeton, USA: Princeton University Press, 1951. 279 p. https://www.jstor.org/stableZj.ctt1b9rzzn

20. Simakov S. S. Modern methods of mathematical modeling of blood flow using reduced order methods. Comput. Res. Model., 2018. Vol. 10, No. 5. P. 581-604. (in Russian) DOI: 10.20537/2076-7633-2018-10-5-581-604

21. Solonnikov V. A. On the solvability of boundary and initial-boundary value problems for the Navier-Stokes system in domains with noncompact boundaries. Pacific J. Math., 1981. Vol. 93, No. 2. P. 443-458. URL: https://projecteuclid.org/euclid.pjm/1102736272

22. Van Dyke M. Perturbation Methods in Fluid Mechanics. New York, London: Academic Press, 1964. 229 p.