A Note on Kirchhoff’s Junction Rule for Power-Law Fluids

1 Introduction

Non-Newtonian fluid flows through tubes that are either thin or long have many industrial and biological applications. It is mainly due to the fact that numerous real fluids (such as polymer melts and solutions, blood, pulps, cement, and mud) cannot be successfully described by the linear Newton’s law because its viscosity significantly changes with the shear rate (see, e.g. [1]). The simplest way to capture those phenomena is to introduce the shear-dependent viscosity. In this article, we study the case where the viscosity is given by the power law (or Ostwald–de Waele law):

(1)η(u) = μ|D(u)|r − 2. (1)

Here, u denotes the velocity, while D(u)=1/2[∇u + (∇u)^T] is the strain rate tensor, i.e. the symmetric component of the velocity gradient. Two positive constants μ and r in (1) are called the flow consistency and the flow behaviour index, respectively. Depending on the value of the index r, we have shear-thinning (or pseudo-plastic) behaviour for 1<r<2 or shear-thickening (or dilatant) behaviour for r>2. Note that the case r=2 corresponds to the standard Newtonian behaviour. The power-law model has been extensively used in the engineering community because of its simplicity and ability to fit the experimental viscosity curves for many real fluids.

Our goal is to study a multiple pipe system filled with power-law fluid obeying (1). Inspired by numerous applications where several pipes are interconnected forming a branching structure, here we focus at the individual section of the network containing one junction point and m straight pipes with (possibly) different profiles (see Fig. 1). The same arguments of the presented analysis can be employed in the case of a general network with multiple junctions as well. Frequently, different scales naturally appear in these systems, e.g. if the extension in one space direction is much larger than in the other two directions. Thus, we introduce a small parameter ε representing the ratio between the cross-section area and the pipes’ length and find the effective behaviour of the flow via asymptotic analysis with respect to ε.

Figure 1:

Junction of thin pipes (for m=3).

The flow of classical Newtonian fluid through a multiple pipe system has so far been addressed in several papers; we refer the reader to [2–5]. The case of power-law flow through a single pipe has been extensively studied for various regimes of the flow, see, e.g. [6–9]. In particular, starting from the Stokes system with a shear-dependent viscosity, the power-law Poiseuille flow (PLPF) has been rigorously derived in [10] via two-scale convergence. Here, we formally confirm that in each pipe of our multiple pipe system, the flow remains in the form of PLPF, at least far from the junction (see Section 3). However, to determine the complete form of PLPF, one needs to deduce the pressure drop along the pipe, i.e. the value of the pressure at the junction point.¹ Using the incompressibility of the flow, we obtain the equation satisfied by the junction pressure p₀, see Section 4. In analogy with electrical circuits, we call it Kirchhoff’s junction rule or simply Kirchhoff’s law. In a case of a simple system consisting of two pipe connected at a junction, it turns out that Kirchhoff’s law can be explicitly solved yielding the explicit formula for p₀. In general situation (the system of m pipes, m>2), we cannot derive the explicit solution of the obtained equation due to its nonlinearity. Nevertheless, we can prove that it generates the unique solution p₀ and that represents our main contribution. Moreover, in the vicinity of junction, the interior layer has been observed in which the fluid exhibits behaviour different than far from the junction. For that reason, in Section 5, we correct the asymptotic approximation for dilatant flows by introducing non-Newtonian Leray’s problem and performing the appropriate matching procedure. Last but not least, the error analysis is provided in Section 6, indicating the order of accuracy of our method. To conclude, we believe that the presented result could be instrumental for creating more efficient numerical algorithms acknowledging the non-Newtonian effects on the pipe networks flow.

2 Description of the Problem

Let Sεi ⊂ R2 be a smooth-bounded domain denoting the cross section of the i^th pipe, i=1, …, m, m ∈ N. The i^th pipe of the system with length ℓ_i in (possibly) different coordinate systems (eji)j = 1, 2, 3 is then defined as

Ωεi = {x = (x1i, x∗i)∈R3 : 0 < x1i < ℓi, x∗i = (x2i, x3i)∈Sεi},i = 1, … m.

Motivated by the applications, we want to study the system in which the pipes are either very thin or very long. The usual way to describe that is to introduce the small (dimensionless) parameter ε as the ratio between pipe’s diameter and its length and put Sεi = εSi. Our thin (or long) pipes are connected at junction point O where we choose to put the origin. In view of that, we introduce Ωε0 = εΩ0 as a small reservoir where the pipes meet, with Ω⁰ being a smooth-bounded domain containing O. The domain under consideration is then defined as (see Fig. 1)

(2)Ωε = ∪i = 0mΩεi. (2)

The end of the i^th pipe is denoted by

Σεi = {(ℓi, x∗i) ∈ R3 : x∗i = (x2i, x3i) ∈ εSi},   i = 1,…,m.

Finally, the lateral boundary of Ω_ε is denoted by Γε = ∂Ωε \ (∪i = 1mΣεi).

We suppose that the pipe system Ω_ε is filled with non-Newtonian fluid with viscosity η obeying the power law:

η(u) = μ|D(u)|r − 2,  1 < r < + ∞ ,  μ = const.  > 0.

Here, D(u)=1/2[∇u + (∇u)^T] is the symmetrised gradient of the velocity, namely

D(u)jk = 12(∂uj∂xk + ∂uk∂xj),  |D(u)| = (∑j, k = 13|D(u)jk|2)1/2.

The equations to be considered have the following form:

(3)− μ   div {|D(uε)|r − 2D(uε)} + ∇pε = 0 in  Ωε, (3)

(4)div  uε = 0  in  Ωε. (4)

Remark 1As you can see, we assume that the flow is sufficiently slow so that the inertial term can be neglected. Such assumption is needed to have the well-posed problem²and not for the sake of our analysis. Indeed, if we write (3) in non-dimensional form and employ the generalised Reynolds number Re_gderived for cylindrical ducts of arbitrary cross section (see [11]), we get Re_g= O (ε^γ) with γ being such that the inertial term does not contribute to the macroscopic model.

To complete the problem, suitable boundary conditions have to be added. We impose a standard no-slip condition for the fluid velocity on the lateral boundary and prescribe the values of pressures p_i at the end of each pipe. In view of that, the boundary conditions read as follows:

(5)uε = 0  on  Γε, (5)

(6)uε  ×  e1i = 0 ,   pε = pi  on  Σεi,  i = 1,…,m. (6)

Our goal is to investigate the effective behaviour of the flow described by (3–6), as ε → 0.

3 Power-Law Poiseuille Flow

Let r′ be the conjugate exponent of r, i.e. r′=r/(r – 1). It is reasonable to expect that in each pipe Ωεi of our system, the flow remains in the so-called PLPF, at least far from the junction. Indeed, by substituting

(7)uε = εr′Ui(x1i, x∗iε),  pε = Pi(x1i, x∗iε) + εQi(x1i, x∗iε) (7)

in the governing equations, we arrive at

(8)ε− 1 :    ∇y∗iPi = 0   ⇒   Pi = Pi(x1i), (8)

(9)1 :    − μ   div y∗i {|Dy∗i(Ui)|r − 2Dy∗i(Ui)} + ∇y∗iQi = −  ∂Pi∂x1i e1i, (9)

(10)εr′ − 1 :  div y∗i Ui = 0 . (10)

The scaling in (7) is suggested by the result from [10] and will be rigorously justified in Section 6 by proving the error estimates. Here, we introduce the fast variable y∗i = x∗i/ε, so the above equations are posed in ε-independent domains Ωi = {(x1i, y∗i) ∈ R3  :  0 < x1i < ℓi,  y∗i = (y2i, y3i) ∈ Si}. We also employ the following notations for the formal partial differential operators:

∇y∗iΦi = ∂Φi∂y2i e2i + ∂Φi∂y3i e3i ,   div y∗iVi = ∂V2i∂y2i + ∂V3i∂y3i ,  Vji = Vi ⋅ eji, j = 1, 2, 3,Dy∗i(Vi)jk = 12(∂Vji∂yki + ∂Vki∂yji), j, k = 2, 3, Dy∗i(Vi)1j = Dy*i(Vi)j1 = 0, j = 1, 2, 3.

Taking into account the no-slip boundary condition for the velocity and the values of the prescribed pressures at pipes ends, the system (8–10) can be solved by taking

Qⁱ=0 and

(11)Ui(y∗i) = wi(y∗i)|pi − p0ℓi|r′ − 2(p0 − piℓi)e1i, y∗i∈Si, (11)

(12)Pi(x1i) = p0 + (pi − p0ℓi)x1i,  x1i∈(0, ℓi). (12)

Here, p₀ denotes the (constant) value of the pressure at the junction point O, while wⁱ is the solution of the auxiliary problem posed on the cross section of the i^th pipe, namely:

(13)− μ div  (|∇ωi|r − 2∇ωi) = 2r/2 in Si, ωi=0 on  ∂Si. (13)

By a simple integration by parts, we deduce

(14)〈ωi〉 :  = ∫Siωi dy∗i = μ2r/2∫Si|∇ωi|r dy∗i > 0. (14)

Remark 2If the pipes are assumed to have circular cross sections, namely Sⁱ= {y ∈ R²: |y|<R_i} (R_i=const.>0), then wⁱcan be explicitly computed as

(15)ωi(y∗i) = 12r/2 − 1μr − 1r′(Rir′ − |y∗i|r′). (15)

Then, in view of (11), the non-zero component of the velocity in the i^thpipe has the form

Uxi(y∗i) = 12r′/2 − 1μr′ − 1r′(Rir′ − |y∗i|r′)|pi − p0ℓi|r′ − 2(p0 − piℓi),y∗i∈Si.

The obtained result is most conveniently illustrated using the dimensionless quantities defined by

U^xi = UxiVmaxi,   y^∗i = y∗iRi,

whereVmaxiis the maximum velocity (obtained fory2i = y3i = 0). The dimensionless velocity fieldU^xi = 1 − |y^∗i|r/(r − 1)is plotted in Figure 2 for different values of the flow behaviour index r. For r=2 (red line) we observe the standard parabolic profile of Newtonian fluid. For r<2 (i.e. for pseudo-plastic fluids), we see that the velocity profile becomes more flattened in the centre as r decreases. The observed phenomena, namely the flattening in the vicinity of the pipe’s centre line and relegation of the majority of shear toward the pipe’s walls could be indicative for practical implications, for instance, in polymer processing flows (see, e.g. [12]).

Figure 2:

Dimensionless velocity distribution U^xi as a function of |y^∗i| for different values of r.

4 Kirchhoff’s Junction Rule

In order to know the complete form of the PLPF derived in the previous section, note that we still need to determine the value of the junction pressure p₀. Denoting Uεi = εr′Ui(x∗i/ε), we can use the incompressibility of our flow and (11) to obtain

(16)0 = ∑i = 1m∫ΣεiUεi⋅e1i dx∗i = − εr′ + 2∑i = 1m〈ωi〉|pi − p0ℓi|r′ − 2(pi − p0ℓi) (16)

Let us introduce

(17)Hi: = |〈ωi〉|r−2〈ωi〉 > 0 (17)

It can be easily verified that

|Hi|r′ − 2Hi = 〈ωi〉

so that from (16) it follows:

(18)∑i = 1m|Hiℓi(pi − p0)|r′ − 2Hiℓi(pi − p0) = 0. (18)

If we consider a pipe system with only two pipes connected at the junction (i.e. m=2), then we can explicitly compute the junction pressure p₀. Indeed, from (18), we deduce

|H1ℓ1(p1 − p0)|r′ − 2H1ℓ1(p1 − p0) = − |H2ℓ2(p2 − p0)|r′ − 2H2ℓ2(p2 − p0).

The above equality is going to be fulfilled if and only if

H1ℓ1(p1 − p0) = − H2ℓ2(p2 − p0)

implying

p0 = H1ℓ1p1 + H2ℓ2p2H1ℓ1 + H2ℓ2.

Taking into account (17), we get the explicit formula for p₀

(19)p0=∑i = 12piℓi|〈ωi〉|r − 2〈ωi〉∑i = 121ℓi|〈ωi〉|r−2〈ωi〉 (19)

In the case of Newtonian fluid (r=2), (19) reduces to standard Kirchhoff’s law derived in [2].

Remark 3Let us compute the formula (19) in the case of two circular pipes, i.e. for Sⁱ={y ∈ R²: |y|<R_i}, i=1, 2…. Taking into account (15), we get

〈ωi〉 = π2r′/2 − 1μr′ − 1(r′ + 2)Rir′ + 2

leading to

(20)p0 = p1ℓ1R13r − 2 + p2ℓ2R23r − 21ℓ1R13r − 2 + 1ℓ2R23r − 2 (20)

To illustrate the result, we investigate how the value of p₀changes with respect to the flow behaviour index r. In view of that, the junction pressure p₀is plotted in Figure 3 as a function of the index r. A wide range for r is taken (1<r<3), and results for different magnitudes of prescribed pressures p₁and p₂are presented.³We observe an interesting phenomena: if characteristic dimension of one pipe is much larger than the other one (green and yellow lines), the junction pressure significantly changes for 1<r<2 and then stabilises for r>2. It means that, for dilatant flows, the junction pressure remains (almost) the same no matter what r we choose. On the other hand, if we have two pipes of similar dimensions, no such phenomena are observed.

Figure 3:

The junction pressure p₀ as a function of the index r.

In general case (for arbitrary, m>2), the junction pressure p₀ is given by (18). Evidently, it is not reasonable to expect that we can explicitly solve such nonlinear algebraic equation, as we managed in the case m=2. The best we can do is to prove that there exists a unique p₀ such that

(21)F(p1, p2, … ,pm, p0) = 0, F(p1, p2, … ,pm, p0) = ∑i = 1m|Hiℓi(pi − p0)|r′ − 2Hiℓi(pi − p0). (21)

For given pressures p₁, p₂,…, p_m, we consider a function t↦G(t) where

(22)G(t) = F(p1, p2, … ,pm, t). (22)

Note that

(23)limt → + ∞G(t) = −∞, limt → − ∞G(t) = + ∞, (23)

so to prove the existence, it is sufficient to show that t↦G(t) is a strictly decreasing function. To accomplish that, we employ the following inequalities:

(24)for r′ > 2  : C1(|b|r′−2 + |a|r′−2)|b − a|2 ≤ (|b|r′ − 2b − |a|r′ − 2a)⋅(b − a) , (24)

(25)for 1 <r′ < 2  : C2 |b − a|2|b|2 − r′+ |a|2 − r′ ≤ (|b|r′ − 2b − |a|r′ − 2a)⋅(b − a) . (25)

The above inequalities hold true for any a, b ∈ Rⁿ, where C₁ and C₂ denote positive constants independent of a and b. For detailed proof of (24) and (25), we refer the reader to [13], Lemma 3.1. Depending on the value of the exponent r from the power law, we deduce as follows:

• Case 1<r<2 (i.e. r′>2):

  (t − s)[G(t) − G(s)]== − ∑i = 1m{|Hiℓi(pi − t)|r′ − 2Hiℓi(pi − t)−|Hiℓi(pi − s)|r′ − 2Hiℓi(pi − s)}⋅(s − t)≤ − C1∑i = 1m{|Hiℓi(pi − t)|r′ − 2+|Hiℓi(pi − s)|r′ − 2}Hiℓi⋅|s − t|2= − C1∑i = 1m〈ωi〉ℓi{|pi − tℓi|r′ − 2 + |pi − sℓi|r′ − 2}⋅|s − t|2 ≤ 0,for any t, s∈R.

  Obviously,

  (t − s)[G(t) − G(s)] = 0  ⇔ t = s,

  so we conclude that G is a strictly decreasing function. Therefore, due to (23), there exists t ∈ R such that G (t)=0. That proves the existence of the junction pressure p₀ satisfying (18). To prove the uniqueness, we assume the opposite. Then, there exist p₀ and q₀ such that G(p₀)=G(q₀)=0. But then (p₀ – q₀)[G(p₀) – G(q₀)]=0 implying p₀=q₀.

• Case r>2 (i.e. 1<r′<2):

  Instead of (24), here we use the inequality (25) to obtain

  (t − s)[G(t) − G(s)]=≤ − C2∑i = 1m1|Hiℓi(pi − t)|2 − r′ + |Hiℓi(pi − s)|2 − r′Hiℓi⋅|s − t|2= − C2∑i = 1m〈ωi〉ℓi1|pi − tℓi|2 − r′ + |pi − sℓi|2 − r′⋅|s − t|2 ≤ 0, for any t, s∈R.

  Using the same arguments as above, we deduce the existence and uniqueness of p₀ in this case as well.

5 Interior-Layer Correction

From the physical point of view, one should expect that the fluid behaviour is somewhat different in the vicinity of junction than far from it. Thus, in this section, we construct the interior-layer corrector capturing those effects and, consequently, improving the order of accuracy of our approximation.

We start by introducing the interior layer by

Gε = Ωε0∪(∪i = 1mGεi),

where

Gεi = {(x1i, x∗i) ∈ Ωεi : 0 < x1i < − d ε ln ε}.

Parameter d>0 describes the thickness of the interior layer and, as such, it will play an important role later on in the error analysis (see Sec. 6). If we pass to the fast variable y=x/ε, namely if we substitute uε(x) = εr′V(x/ε),p_ε(x)=εΠ(x/ε) into (3) and (4), we obtain

(26)− μ div {|D(V)|r−2D(V)} + ∇Π = 0 in G∞, (26)

(27)div V = 0 in G∞. (27)

The above equations are posed in the domain G∞ = Ω0∪∪i = 1mΩ∞i consisting of a bounded reservoir Ω⁰ and m infinitely long pipes Ω∞i = {y = (y1i, y∗i) ∈ R3 : y1i > 0, y∗i = (y2i, y3i) ∈ Si}. Since far from the reservoir, the flow in pipes remains in the form of PLPF (see Sec. 3), we endow (26) and (27) with the conditions:

(28)V = 0 on ∂G∞, (28)

(29)|y| → + ∞limV = Ui in Ω∞i, i = 1, …, m, (29)

where Uⁱ is given by (11). The Stokes problem (corresponding to the value r=2 in our case) in an infinite domain with asymptotic Poiseuille velocities is called Leray’s problem, and an overview of the related mathematical theory can be found in [14]. Here, the system (26–29) forms its non-Newtonian generalisation and has been theoretically analyzed in [15] by the first author of the present article. In the case of dilatant flow (r>2), the existence and uniqueness of the solution have been established as well as the exponential decay estimate for the solution. To our knowledge, the similar result for pseudo-plastic flows (1<r<2) is still missing and cannot be found in the existing literature. For that reason, we are in position to define the correction in the interior layer G_ε only for r>2. Namely, it has the form

(30)Vε(x) = εr′V(xε), Πε(x) = εΠ(xε), (30)

with (V, Π) being the solution of the non-Newtonian Leray’s problem (26–29).

Note we still need to perform some matching with the outer flow. Indeed, the PLPF (Uεi, Pi) representing the solution in

ωεi = {(x1i, x∗i) ∈ Ωεi : − d ε ln ε < x1i < ℓi}

and the interior-layer corrector (Vεi, Πεi) in G_ε are only matched asymptotically (for ε → 0), due to (29). Thus, on each interface,

γεi = {(x1i, x∗i) ∈ Ωεi : x1i = − d ε ln ε}

between G_ε and ωεi, our asymptotic approximation is not continuous but it holds

Vε = Uεi + Eεi,  Eεi(x*i) = εr′[V(− d  ln ε, x∗iε) − Ui(x˜iε)].

Thanks to the exponential decay result of the Leray solution (see [15], theorem 5.1), the residual Eεi can be controlled. Moreover, by a simple change of variables from (29), it follows

∫γεiEεi⋅e1i dx∗i = εr′ + 2∫Si[V(− d ln ε, y∗i) − Ui(y∗i)]d y∗i = 0.

That enables us to introduce hεi as the solution of the standard divergence problem (see, e.g. [14], chapter III.3):

(31)div hεi = 0 in ωεi, hεi = − ε− r′ Eεi on γεi, hεi = 0 on ∂ωεi\γεi. (31)

Finally, we define our (matched) asymptotic solution as

(32)Uεapprox = εr′Uεapprox, where Uεapprox = {Ui(x∗iε)+hεi in ωεiV(xε) in Gε, (32)

(33)Pεapprox = {Pi in ωεiΠε in Gε (33)

6 Error Analysis

We restrict ourselves only on the case r>2. The error estimate can be derived in the similar way as in [2], using, in addition, the monotonicity of the differential operator given by the power law. Indeed, it can be verified by direct computation that, for some C, c>0 (independent on ε) it holds

(34)c∫Ωε|D(Rε)|r ≤ μ ∫Ωε[|D(uε)|r − 2D(uε) − |D(Uεapprox)|r − 2D(Uεapprox)]D(Rε) = = ∑i = 1m∫γεiJεi Rε ≤ C|ε  log  ε|1 − 1/r∑i = 1m|Jεi|Lr′(γεi)|D (Rε)|Lr(Ωε), (34)

where Rε = uε − Uεapprox and Jεi, in i^th pipe, is the jump of the normal component of the stress tensor over the surface γεi. Essentially, we have

Jεi = ε[Π(d|ln  ε|, x∗iε) − d|ln  ε|pi − p0ℓi+ μ|D(V(d|ln  ε|, x∗iε))|r − 2D(V(d|ln  ε|, x∗iε))− μ|D(Ui(x∗iε))|r − 2D(Ui(x∗iε))] ei =  O(εd + 1).

Taking into account that the measure of the surface |γεi| = O(ε2) gives

|D(Rε)|Lr(Ωε) ≤ Cε(4 + d− 3r)/(r−1)|ln ε|1/r.

Thus, following the similar procedure as in [2], we obtain for r′=r/(r – 1):

(35)1|Ωε|∫Ωε|D(uε) − D(Uεapprox|r ≤ Cεβ|ln ε|,β = r(r − 1)(4 + d − r) − 3r − 1− 2. (35)

(36)1|Ωε|∫Ωε|uε − Uεapprox|r ≤ Cεβ+r|ln ε| (36)

(37)1|Ωε|∫Ωε|pε − 𝒫εapprox|r′ ≤ Cεσ|ln ε|1/(r − 1)2, σ = r(4 + d) − 3(r − 1)3 − r′ − 2. (37)

For that velocity estimate to make sense, we need to have β + r>0, which is true for any r>2 and d>0. For the pressure estimate, we need σ>0, leading to the more restrictive condition

d > 1r[(r − 1)2(3r − 2) + 3] − 4.

Thus, the thickness of the interior layer is depending on the flow index r and increases with it. The magnitude of the interior-layer corrector for the velocity is of the order ε^1/r (in L^r norm) and for the pressure ε1+1r′, in L^r′ norm. So, the correctors improve our approximation for the velocity only in the case d>r – 1. On the contrary, for the pressure, the corrector makes the difference only if d>5r² – 7r + 1 – 3/r, which is very large for r>2. Thus, it is reasonable to keep the velocity corrector if we want better approximation, while the pressure corrector is just a technical aid needed for deriving the satisfactory error estimate.

7 Conclusion

The study of non-Newtonian fluid flow has gained much attention in recent years because of its numerous technological applications. Several models of non-Newtonian fluids have been proposed and among all, the simplest and the most common one is the power-law model. On the other hand, engineering practice requires extensive knowledge of the pipe network flow. Therefore, the aim of this article was to consider the steady flow of a power-law fluid through a junction of thin (or long) pipes with possibly different profiles. No serious assumptions are imposed in order to simplify the original 3D problem describing the real physical situation. Using asymptotic analysis with respect to the pipes thickness, the effective behaviour of the flow has been found. In particular, Kirchhoff’s law (18) has been derived yielding the value of the junction pressure. The flow in the vicinity of the junction has been studied in the case of dilatant flow, and the corrected asymptotic approximation (32) and (33) has been built. The accuracy of the obtained approximation has been discussed by providing the corresponding error estimates [see (36) and (37)]. Due to the lack of the existence and decay result concerning the pseudo-plastic flow in unbounded domain, we are not in position to correct the PLPF in the vicinity of the junction for 1<r<2. This is the subject of our current investigation. Nevertheless, we believe that the obtained result is important from both theoretical and practical points of view and, hopefully, could improve the known engineering practice.