Non-Newtonian Flow in Bifurcated Piping and Arterial Networks

Jamie A. Chin , Wilson C. Chin ^*

Stratamagnetic Software, LLC, Houston, Texas, United States

* Correspondence: Wilson C. Chin

Academic Editor: Faik Hamad

Special Issue: Multi-phase Flow with and without Heat Transfer

Received: February 12, 2022 | Accepted: May 06, 2022 | Published: May 23, 2022

Journal of Energy and Power Technology 2022, Volume 4, Issue 2, doi:10.21926/jept.2202018

Recommended citation: Chin JA, Chin WC. Non-Newtonian Flow in Bifurcated Piping and Arterial Networks. Journal of Energy and Power Technology 2022; 4(2): 018; doi:10.21926/jept.2202018.

© 2022 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction to Bifurcated Flows

Engineers have long recognized the important roles played by rheology, pressure, viscous stress and geometry in design. Solutions to bifurcated flows are well known in engineering. For example, in two-branch systems where all three inlet and outlet paths are nearly parallel, as in tuning forks, losses at junctions are small – an elementary solution, offered at the tutorial website www.chegg.com, assumes ideal fluids where total energy is conserved [1]. In process applications, piping networks follow energy laws with empirical “head loss” factors for fittings, roughness and bends. Sometimes velocity profiles (from simpler fluid models) are used to evaluate integrals that estimate total losses. Again, ideal “total energy conservation” and empirical “head loss coefficient” models cannot always be meaningfully used – key issues cannot be understood without analytical tools based on accepted rheological models. For example, interpretation should depend explicitly on measurable parameters related to Newtonian, Power Law, Bingham Plastic or Herschel-Bulkley inputs. More sophisticated three-dimensional finite element methods write algebraic equations to large numbers of nodal points that are then solved. Despite the apparent generality, their time-consuming solutions do depend on meshing details – and while predictive capabilities are possible, they will require careful calibration and specific empirical inputs. In all of these methods, physical dependencies on fluid and geometric parameters can be hidden or lost, thus limiting their use as diagnostic or design tools. Simple questions related to the roles of radius, “n” and “k” have not been be straightforwardly answered.

Bifurcation applications are broad and not restricted to engineering design. Medical research, increasingly sophisticated, has challenged physicists and rheologists in recent years. Importantly, medical researchers and practitioners study health problems using the same parameters – the rheological models cited earlier are able to describe a broad range of multiphase effects inherent in heterogeneous blood flows. For instance, it is now recognized that the main arteries host largely Newtonian flows, while complicated non-Newtonian effects are relegated to important secondary systems with smaller diameters. At the finest capillary level, nonlinear effects are even more pronounced. These massive bifurcated systems often defy analysis and do not allow detailed description. Fortunately, the principal flows branching from the main aorta that control overall bodily functions can be modeled – these branched configurations also form the focus of this investigation.

Models that yield physical relationships between volume flow rate, pressure, viscous stress and geometric parameters would promote more optimized pipe flows as well as an improved understanding of degenerative health conditions. In this paper, we provide simple yet new and useful solutions to important practical bifurcation problems appearing in both engineering practice and medical research. The analytical models and numerical implementation, as well as their extensions to handle clogged flows with arbitrary cross-sectional geometries, are described in detail with Fortran source code offered where possible. The role of viscous shear stress in dislodging clogs and in preventing their formation cannot be de-emphasized. Hopefully, the approach to the overall problem presented here will prove immediately useful to practitioners in the field. Additionally, simple interactive software screens that host the new models are shown that support ease-of-use and rapid analyses. We refer the reader to Chhabra and Richardson [2] and Schlichting [3] for basic fluid dynamic models, and to Hussain, Kar and Puniyani [4], Lee, Xue, Nam, Lim and Shin [5], Sochi [6] and Yamamoto et al. [7] for biological applications.

1.1 Bifurcation Model and Analytical Approach

Consider a system with a single pipe (supporting an inlet flow with given flow rate) that “bifurcates” or branches into N-1 multiple conduits (hosting N-1 outlet flows with prescribed exit pressures). It is physically clear, provided each piping element is longer than the dimensions of the junction, that secondary and turbulent flow effects will be localized at the intersection. Thus, the analytical flow rate formula for each constituent pipe flow would take a form identical to that of the “single pipe alone” problem for that particular pipe.

Only the inlet and outlet pressures for each of the N conduits need to be specified. All of these pressures would be available if the junction pressure can be determined. We obtain this value by assuming total mass conservation – this implies that total energy and momentum are not conserved, which is anticipated in real world applications (junction losses, assumed to be small, are excluded here). Hence, relationships connecting flow rate, axial velocity, pressure, apparent viscosity, shear stress and shear rate are easily calculated from the underlying rheological solution.

By contrast, methods taught in elementary fluid mechanics encompass a broad range of approximations. The simplest ideal flow models assume unrealistic total energy conservation – this unfortunately suggests that power is not required to maintain the flow [1]. On the other hand, the empirical “head loss” methods used in the piping industry only utilize measured friction data dependent on surface roughness, inlet conditions, fittings details, local bends and so on. Rheology is often not considered. These methods are unlikely to produce useful predictive results in more general flow problems.

Our use of local “single pipe solutions” for constituent pipes in steady bifurcated systems is somewhat analogous to the WKB methods developed in transient wave propagation [8]. There, large space-time solutions for complicated propagation fields are assumed in a form “ a sin (kx - ωt)” similar to exact elementary solutions, except that amplitude, wavenumber and frequency are now slowly varying functions of space and time. It is important to understand that our use of single pipe solutions is not restricted to those having circular cross-sections. Our general approach applies to square, triangular, highly clogged or even mixed combinations of these (for any host rheology) provided an analytical representation for each of the constituent piping elements is available. Moreover, the method is applicable to more complicated situations where body forces, such as gravitational and electrical, are found.

1.2 Different Rheological Applications

Our method applies to different types of liquid stress-strain constitutive relations. In the work presented here, we illustrate four popular models, namely Newtonian, Power Law, Bingham Plastic and Herschel-Bulkley fluids. The four “single pipe alone” circular tube solutions used to host bifurcation models in each case are listed below where the constants are defined separately in this write-up. These models, in the order listed, are shown to require added numerical support as physical complexity increases.

• Newtonian, $\mathrm{Q}=\unicode[Times]{x03C0} \mathrm{R}^{4} \Delta \mathrm{p} /(8 \unicode[Times]{x03BC} \mathrm{L})$
• Power Law, $\mathrm{Q} /\left(\unicode[Times]{x03C0} \mathrm{R}^{3}\right)=[\mathrm{R} \Delta \mathrm{p} /(2 \mathrm{~kL})]^{1 / \mathrm{n}} \mathrm{n} /(3 \mathrm{n}+1)$
• Bingham Plastic, $\mathrm{Q} /\left(\unicode[Times]{x03C0} \mathrm{R}^{3}\right)=\left[\unicode[Times]{x03C4}_{\mathrm{w}} /(4 \unicode[Times]{x03BC})\right]\left[1-4 / 3\left(\unicode[Times]{x03C4}_{\mathrm{O}} / \unicode[Times]{x03C4}_{\mathrm{w}}\right)+1 / 3\left(\unicode[Times]{x03C4}_{\mathrm{O}} / \unicode[Times]{x03C4}_{\mathrm{w}}\right)^{4}\right]$
• Herschel-Bulkley, $\mathrm{Q} /\left(\unicode[Times]{x03C0} \mathrm{R}^{3}\right)=\mathrm{k}^{-1 / \mathrm{n}}(\mathrm{R} \Delta \mathrm{p} / 2 \mathrm{L})^{-3}\left(\mathrm{R} \Delta \mathrm{p} / 2 \mathrm{L}-\unicode[Times]{x03C4}_{\mathrm{O}}\right)^{(\mathrm{n}+1) / \mathrm{n}} \times$
  $\left[\left(\mathrm{R} \Delta \mathrm{p} / 2 \mathrm{L}-\unicode[Times]{x03C4}_{\mathrm{O}}\right)^{2} \mathrm{n} /(3 \mathrm{n}+1)+2 \unicode[Times]{x03C4}_{\mathrm{O}}\left(\mathrm{R} \Delta \mathrm{p} / 2 \mathrm{L}-\unicode[Times]{x03C4}_{\mathrm{O}}\right) \mathrm{n} /(2 \mathrm{n}+1)+\unicode[Times]{x03C4}_{\mathrm{O}}^{2} \mathrm{n} /(\mathrm{n}+1)\right]$

These choices are directed at highlighting custom techniques specific to different fluid models. In the simplest linear Newtonian case, it is possible to obtain closed form, analytical solutions because of algebraic simplicities. For Power Law flows, the nonlinearity disallows elegant solutions, and a “half-step” iteration is instead required to solve for junction pressure. For Bingham Plastics and Herschel-Bulkley fluids with yield stress, further complications arise because expressions for the inlet pressure cannot be obtained. For these two problems, two successive half-step iteration procedures are needed to fulfill our solution objectives. Other possible rheologies not considered here, for example, include

• Ellis, $\unicode[Times]{x03C4}=-\mathrm{d u} / \mathrm{d r} /\left(\mathrm{A+B} \unicode[Times]{x03C4}^{\alpha-1}\right)$
  $\mathrm{Q} / \unicode[Times]{x03C0} \mathrm{R}^{3}=\mathrm{A} \unicode[Times]{x03C4}_{w} / 4+\mathrm{B} \unicode[Times]{x03C4}_{w}{ }^{\alpha} /(\alpha+3)$
              $=\mathrm{A}(\mathrm{R} \Delta \mathrm{p} / 2 \mathrm{L}) / 4+\mathrm{B}(\mathrm{R} \Delta \mathrm{p} / 2 \mathrm{L})^{\alpha} /(\alpha+3)$
• Casson, with constants $\unicode[Times]{x03C4}_{\mathrm{O}}, \unicode[Times]{x03BC}_{\mathrm{O}}$ and a dimensionless  $\xi=\unicode[Times]{x03C4}_{\mathrm{O}} /\left\{\frac{\left(\mathrm{P}_{\mathrm{O}}-\mathrm{P}_{\mathrm{L}}\right) \mathrm{R}}{2 \mathrm{L} .}\right\}, $
  $\mathrm{Q}=\unicode[Times]{x03C0} \mathrm{R}^{4}\left(\mathrm{P}_{\mathrm{O}}-\mathrm{P}_{\mathrm{L}}\right)\left\{1-(16 / 7) \xi^{1 / 2}+(4 / 3) \xi-(1 / 21) \xi^{4}\right\} /\left(8 \unicode[Times]{x03BC}_{\mathrm{O}} \mathrm{L}\right)$
• Eyring, with constants $\unicode[Times]{x03C4}_{\mathrm{O}}, \lambda_{O}$ and $\unicode[Times]{x03C4}_{\mathrm{R}}=\left(\mathrm{P}_{\mathrm{O}}-\mathrm{P}_{\mathrm{L}}\right) \mathrm{R} / 2 \mathrm{L}$
  $\mathrm{Q}=2 \unicode[Times]{x03C0} \mathrm{R}^{3}\left(\frac{\unicode[Times]{x03C4}_{_{O}}}{\unicode[Times]{x03C4}_{_{R}}}\right)^{3}\left[\left\{1 / 2\left(\unicode[Times]{x03C4}_{R} / \unicode[Times]{x03C4}_{\mathrm{O}}\right)^{2}+1\right\} \cosh \left(\unicode[Times]{x03C4}_{R} / \unicode[Times]{x03C4}_{\mathrm{O}}\right)-\left(\unicode[Times]{x03C4}_{R} / \unicode[Times]{x03C4}_{\mathrm{O}}\right) \sinh \left(\unicode[Times]{x03C4}_{R} / \unicode[Times]{x03C4}_{\mathrm{O}}\right)-1\right] / \lambda_{O}$  

1.3 Validation Procedures

Detailed discussions for the above formulas and for other rheologies are available in excellent references, for example, Bird, Armstrong and Hassager [9], Bird, Stewart and Lightfoot [10], Ferry [11] and other classic rheology books. We will not rederive (except use) the final results. We do emphasize several editorial decisions made in presenting our bifurcation models.

• We are not providing exhaustive parameter studies for particular models, since these would dilute the impact of our focus on methodology. For example, even in the simple Power Law limit, variations in the parameter ranges of (n, k) are broad and depend on the industry and application (as such, n < 1 suffices in fulfilling our objectives and n > 1 is unnecessary).
• It has been suggested that experimental validation is required. This is not so, since extraneous effects related to inlet conditions, surface roughness, junction loss and more will not be consistent with the model assumptions applied. Assessments for different models must be obtained using methods that rely solely on the math assumptions taken, as explained next.
• Our validations are comprehensive. (1) We start with our analytical Newtonian bifurcation solution, verify expected physical dependencies related to pressure drop and viscosity, and ensure that flow symmetries are obtained for symmetrical bifurcation placements. (2) For the Power Law solution, we operate the iterative algorithm for values of the dimensionless n that are very close to unity with k=µ (the Newtonian viscosity) to ensure agreement with the Newtonian solution in (1). (3) For Bingham Plastics with n = 1 and nonzero yield stress, the “two half-step” iteration is operated with almost vanishing yield stress to verify that results agree with the “single half-step” approach for Power Law fluids. (4) Finally, for Herschel-Bulkley fluids, the yield stress code is operated to ensure that results are consistent with the Bingham Plastic model in (3) with n = 1. 
• In summary, the objective of the work presented here focuses on the mechanics of the solution strategy and its implementation for different rheologies. Thus, tables of numerical results are not offered as they will be very limited in usefulness. Instead, the algorithms and user interfaces themselves are described in detail. We emphasize that our work relates to both piping and arterial networks in living tissue as both applications are important in seemingly different but actually related disciplines. 

2. Newtonian Fluids in General Circular and Clogged Conduits

2.1 Newtonian Flow in Simple Bifurcations

We will present some useful analytical models for piping systems and blood vessel systems that support Newtonian flows. We first consider fluid flowing through a main conduit that bifurcates into two different branches. Physical assumptions and their mathematical expression are explained at each step of the derivations. These require no more than elementary algebra, but the problem and results do not appear to have been given previously. Practical applications are emphasized in the discussions. Then in subsequent sections, more complicated flow networks are considered, now requiring much less explanation. In many cases, user-friendly software is available to perform calculations rapidly. Extensions to conduits with arbitrary non-circular cross-sections will also be considered – this methodology is relevant to flows in irregular clogged pipes. For instance, calculated viscous stresses at complicated debris-fluid interfaces may support cleaning and remediation efforts.

We start with the Hagen-Poiseuille formula for the volume flow rate “$ \text{Q}=\unicode[Times]{x03C0} \text{ R}^{4}\left(\text{P}_{i} - \text{P}_{o}\right) /(8 \unicode[Times]{x03BC} \text{L})$ .” Its assumptions are numerous and well known: (a) a rigid circular cross-section of radius R, (b) a straight conduit with length L satisfying L >> R, (c) a constant fluid viscosity of μ, (d) end conditions with P_i and P_o respectively being inlet and outlet pressures, (e) steady-state laminar liquid flow, (f) smooth walls and (g) a location away from inlets. Our convention assumes that the volume flow rate is Q > 0 if P_i >P_o . Fluid density does not appear in Hagen-Poiseuille’s law because steady-state conditions are assumed – density, however, is important in transient applications. The two-branch bifurcation shown in Figure 1 applies to both engineering and medical applications.

Figure 1 Two-branch bifurcation (top), multi-branch geometry (bottom).

It is possible to develop extended applications of the classic Hagen-Poiseuille rule. Consider the uneven bifurcation in Figure 1, where the main flow rate Q₁ in “1” feeds into branches “2” and “3” as shown at the top of Figure 1. For a constant density flow, mass conservation requires that Q₁ = Q₂ + Q₃. (Here mass considerations disallow momentum and energy conservation, both of which are not possible in the presence of irreversible viscous wall losses). Applying the rule to each flow segment directly yields $\text{Q}_{1}=\unicode[Times]{x03C0} {\text{R}_{2}}^{4}\left(\text{P}_ \mathcal a-\text{P}_{o, 2}\right) /\left(8 \unicode[Times]{x03BC} \text{L}_{2}\right)+\unicode[Times]{x03C0} \text{R}_{3}{ }^{4}\left(\text{P}_ \mathcal a-\text{P}_{o, 3}\right) /\left(8 \unicode[Times]{x03BC} \text{L}_{3}\right)$. Here, P_o,2 and P_o,3 are branch “2” and “3”’s outlet pressures, while at the junction “$a$,” we determine the pressure as

\[ \mathrm{P}_ \mathcal a=\left\{8 \unicode[Times]{x03BC} \mathrm{Q}_{1} / \unicode[Times]{x03C0}+\left(\mathrm{R}_{2}{ }^{4} \mathrm{P}_{\mathrm{o}, 2} / \mathrm{L}_{2}+\mathrm{R}_{3}{ }^{4} \mathrm{P}_{\mathrm{o}, 3} / \mathrm{L}_{3}\right)\right\} /\left\{\mathrm{R}_{2}{ }^{4} / \mathrm{L}_{2}+\mathrm{R}_{3}{ }^{4} / \mathrm{L}_{3}\right\} \tag{2.1} \]

This can be substituted into the formulas

\[ \text{Q}_{2}=\unicode[Times]{x03C0} {\text{R}_{2}}^{4}\left(\text{P}_ \mathcal a-\text{P}_{o,2}\right) /\left(8 \unicode[Times]{x03BC} \text{L}_{2}\right) \tag{2.2} \]

\[ \text{Q}_{3}=\unicode[Times]{x03C0} {\text{R}_{3}}^{4}\left(\text{P}_ \mathcal a-\text{P}_{o,3}\right) /\left(8 \unicode[Times]{x03BC} \text{L}_{3}\right) \tag{2.3} \]

to produce the volume flow rates Q₂ and Q₃ into branches “2” and “3” respectively.

We could also have expressed “Q₁ = Q₂ + Q₃” in the form R₁⁴ (P_i – P_a) /L₁ = R₂⁴ (P_a – P_o,2) /L₂ + R₃⁴ ( – P_o,3) /L₃ to show P_α= (R₁⁴ P_i/L₁ + R₂⁴ P_o,2 /L₂ + R₃⁴ P_o,3 /L₃)/(R₁⁴ /L₁ + R₂⁴ /L₂ + R₃⁴ /L₃). If we substitute this into $\text{Q}_{1}=\unicode[Times]{x03C0} \text{R}_{1}{ }^{4}\left(\text{P}_{i}-\text{P}_ \mathcal a\right) /\left(8 \unicode[Times]{x03BC} \text{L}_{1}\right)$, we obtain

\[ \mathrm{Q}_{1}=\left\{\unicode[Times]{x03C0} \frac{\mathrm{R}_{1}{ }^{4}}{8 \unicode[Times]{x03BC} \mathrm{L}_{1}}\right\} \times\left[\mathrm{P}_{\mathrm{i}}-\left\{\frac{\frac{\mathrm{R}_{1}{ }^{4} \mathrm{P}_{\mathrm{i}}}{\mathrm{L}_{1}}+\frac{\mathrm{R}_{2}{ }^{4} \mathrm{P}_{o,2}}{\mathrm{L}_{2}}+\frac{\mathrm{R}_{3}{ }^{4} \mathrm{P}_{\mathrm{o}, 3}}{\mathrm{L}_{3}}}{\frac{\mathrm{R}_{1}{ }^{4}}{\mathrm{L}_{1}}+\frac{\mathrm{R}_{2}{ }^{4}}{\mathrm{L}_{2}}+\frac{\mathrm{R}_{3}{ }^{4}}{\mathrm{L}_{3}}}\right\}\right] \tag{2.4} \]

which can be inverted to yield

\[ \mathrm{P}_{\mathrm{i}}= \frac{\mathrm{Q}_{1}+\left\{\unicode[Times]{x03C0} \frac{\mathrm{R}_{1}{ }^{4}}{8 \unicode[Times]{x03BC} \mathrm{L}_{1}}\right\} \left\{\frac{\frac{\mathrm{R}_{2}{ }^{4} \mathrm{P}_{\mathrm{o}, 2}}{\mathrm{L}_{2}}+\frac{\mathrm{R}_{3}{ }^{4} \mathrm{P}_{\mathrm{o}, 3}}{\mathrm{L}_{3}} } {\frac{\mathrm{R}_{1}{ }^{4}}{\mathrm{L}_{1}}+\frac{\mathrm{R}_{2}{ }^{4}}{\mathrm{L}_{2}}+\frac{\mathrm{R}_{3}{ }^{4}}{\mathrm{L}_{3}}}\right\}} {\left\{\unicode[Times]{x03C0} \frac{\mathrm{R}_{1}{ }^{4}}{8 \unicode[Times]{x03BC} \mathrm{L}_{1}}\right\} -\frac{\left\{\unicode[Times]{x03C0} \frac{\mathrm{R}_{1}{ }^{8}} {8 \unicode[Times]{x03BC} \mathrm{L}_{1}{ }^{2}}\right\} } {\frac{\mathrm{R}_{1}{ }^{4}} {\mathrm{L}_{1}}+\frac{\mathrm{R}_{2}{ }^{4}}{\mathrm{L}_{2}}+\frac{\mathrm{R}_{3}{ }^{4}}{\mathrm{L}_{3}}}} \tag{2.5} \]

Finally, the viscous shear stresses $\unicode[Times]{x03C4}$ acting at the walls of branches “1,” “2” and “3” are given by

\[ \unicode[Times]{x03C4}_1=\left(\text{P}_{i}-\text{P}_ \mathcal a\right) \frac{\text{R}_{1}}{2 \text{L}_{1}} \tag{2.6} \]

\[ \unicode[Times]{x03C4}_2=\left(\text{P}_ \mathcal a-\text{P}_{o,2}\right) \frac{\text{R}_{2}}{2 \text{L}_{2}} \tag{2.7} \]

\[ \unicode[Times]{x03C4}_3=\left(\text{P}_ \mathcal a-\text{P}_{o,3}\right) \frac{\text{R}_{3}}{2 \text{L}_{3}} \tag{2.8} \]

which follow from global momentum balances. These equations all solve different flow problems.

2.1.1 Case 1 Flow Rate Q₁ Prescribed

If the main flow rate Q₁ is given, together with the outlet pressures P_o,2 and P_o,3, and if all lengths, radii, and viscosity are specified, then the volume flow rates in bifurcations “2” and “3” are determined by Equations 2.2 and 2.3 where the junction pressure P_a is defined by Equation 2.1. The inlet pressure where the main flow enters is calculated from Equation 2.5. This formulation addresses the following problems. If the pumping is achieved, say with a mechanical device, will the inlet pressure P_i exceed a pre-defined “burst pressure” and lead to wall breakage? Will the rates Q₂ and Q₃ be rapid enough to deliver needed oxygen and nutrients in the case of blood flow? Will viscous shears be so excessive so as to induce long-term wear and tear? How might R and L be altered to achieve safety targets? In health care, options may include increases in R (through clog removal and stent usage), decreases in L (as in a shortening of varicose veins), and reductions in viscosity using blood thinners.

2.1.2 Case 2 Inlet Pressure P_i Prescribed

If the inlet pressure is given, together with the outlet pressures P_o,2 and P_o,3, and if all lengths, radii and viscosity are specified, then the main flow rate Q₁ is known from Equation 2.4. This can be substituted in Equation 2.1 to give the junction pressure P_a, which is in turn used in Equations 2.2 and 2.3 to provide Q₂ and Q₃. We might similarly ask, “Is Q₁ rapid enough to deliver needed oxygen and nutrients?” For the geometries assumed, will Q₂ and Q₃ meet flow requirements? How can the junction pressure P_a be increased to improve these volume flow rates? Can we alter R₂ and L₂ to effect changes in branch “3,” if this branch itself is not easily accessible?

2.1.3 Case 3 Identical Outlet Pressures P_o,2 and P_o,3 Given

In this case, we write P_o,2 = P_o,3 = P_o. Here, the geometries for branches “2” and “3” may differ. For instance, one conduit may be longer than the other, with both having different radii – however, we assume that their outlets are sufficiently close in proximity that their exit pressures will also be close. In other applications, the two may be geometrically identical but situated far apart – an obvious human body example is found in the left and right iliac arteries that transport blood through the pelvis and legs. Then, Equation 2.1 simplifies considerably and substitution into Equations 2.2 and 2.3 shows that

\[ \text{P}_ \mathcal a=\text{P}_{o}+\frac{\left(8 \unicode[Times]{x03BC} \frac{\text{Q}_{1}}{\unicode[Times]{x03C0}}\right)}{\frac{{\text{R}_{2}}^{4}}{\text{L}_{2}}+\frac{{\text{R}_{3}}^{4}}{\text{L}_{3}}} \tag{2.9} \]

\[ \text{Q}_{2}=\frac{\frac{\text{Q}_{1} \text{R}_{2}{ }^{4}}{\text{L}_{2}}}{\frac{\text{R}_{2}{ }^{4}}{\text{L}_{2}}+\frac{\text{R}_{3}{ }^{4}}{\text{L}_{3}}} \tag{2.10} \]

\[ \mathrm{Q}_{3}=\frac{\frac{\mathrm{Q}_{1} \mathrm{R}_{3}{ }^{4}}{\mathrm{L}_{3}}}{\frac{\mathrm{R}_{2}{ }^{4}}{\mathrm{L}_{2}}+\frac{\mathrm{R}_{3}{ }^{4}}{\mathrm{L}_{3}}} \tag{2.11} \]

These results are interesting. First, Q₂ and Q₁ are related only by a simple geometric factor (R₂⁴/L₂)/(R₂⁴/L₂ + R₃⁴ /L₃) that is independent of viscosity and pressure – this similarly applies to Q₃ and Q₁. Second, the quantity “Q₂ + Q₃” can be evaluated from Equations 2.10 and 2.11 to show that Q₂ + Q₃ = Q₁, which states that mass is conserved as assumed.

And third, a simple “energy like” law follows from Equations 2.10 and 2.11. Squaring these leads to Q₂² = (Q₁²R₂⁸/L₂²)/(R₂⁴ /L₂ + R₃⁴ /L₃)² and Q₃² = (Q₁²R₃⁸/L₃²)/(R₂⁴ /L₂ + R₃⁴ /L₃)². Now, adding the two gives

\[ \mathrm{{Q}_{2}}^{2}+\mathrm{{Q}_{3}}^{2}=\mathrm{{Q}_{1}}^{2} \frac{\frac{\mathrm{{R}_{2}}^{8}}{\mathrm{{L}_{2}}^{2}}+\frac{\mathrm{{R}_{3}}^{8}}{\mathrm{{L}_{3}}^{2}}}{\left(\frac{\mathrm{{R}_{2}}^{4}}{\mathrm{{L}_{2}}}+\frac{\mathrm{{R}_{3}}^{4}}{\mathrm{{L}_{3}}}\right)^{2}}=\mathrm{{Q}_{1}}^{2}\left\{\frac{\frac{\mathrm{{R}_{2}}^{8}}{\mathrm{{L}_{2}}^{2}}+\frac{\mathrm{{R}_{3}}^{8}}{\mathrm{{L}_{3}}^{2}}}{\frac{\mathrm{{R}_{2}}^{8}}{\mathrm{{L}_{2}}^{2}}+2 \mathrm{{R}_{2}}^{4} \frac{\mathrm{{R}_{3}}^{4}}{\mathrm{{L}_{2}} \mathrm{{L}_{3}}}+\frac{\mathrm{{R}_{3}}^{8}}{\mathrm{{L}_{3}}^{2}}}\right\} \tag{2.12a} \]

One can think of this as an equation for energy that is dissipated in the bifurcation process, with { } representing the exact loss factor. And since { } is a fraction less than unity, with its denominator exceeding the value of its numerator, we clearly have the inequality

\[ \mathrm{{Q}_{2}}^{2}+\mathrm{{Q}_{3}}^{2}< \mathrm{{Q}_{1}}^{2} \tag{2.12b} \]

Thus, while the mass conservation in Q₂ + Q₃ = Q₁ holds, the sum of the squares in volume flow rates, a quantity related to kinetic energy, is not conserved. Instead, it must follow Equation 2.12b for the Newtonian model assumed. As given, the above relations are useful for quick “back of the envelope” estimates. In closing, we emphasize that the matching condition Q₂ + Q₃ = Q₁ used allows an abrupt lossy transition from a single conduit to a doubly-bifurcated system. In reality, fluid-dynamical motions can be very complicated at such junctions, but these will require lengthy and detailed three-dimensional finite element simulations.

2.2 Software – CODE-1 for Two Uneven Bifurcated Arteries with Q₁ Specified

The results derived above are easily programmed on spreadsheet utilities, but they are more conveniently accessed in the CODE-1 utility of Figure 2. A symmetric artery system is shown, e.g., as for the left and right iliac arteries mentioned earlier, but we emphasize that our model also applies to systems without symmetry (this “app” is also applicable to engineering flows).

Figure 2 Forward model with Q₁ given (Reference, CODE-1).

2.2.1 A Sample Calculation

To validate our bifurcated flow program, we assume that branches “1,” “2” and “3” are geometrically identical. In Figure 2, white text boxes are reserved for data inputs, while gray boxes are intended for displays of non-editable calculated results. Clicking “Find” in Figure 3 automatically computes all of the flow attributes derived above. The program correctly shows that the assumed flow rate of “100 cc/s” splits into two identical flows of 50 cc/s each. The pressure drops in branches “2” and “3” are 2.00373 – 2 = 0.00373 psi, while that in “1” is 2.01119 – 2.00373 = 0.00746 psi. This is exactly twice “0.00373 psi” and is consistent with a flow rate ratio of two, obtained on comparing 100 to 50 cc/s. This linear dependence of volume flow rate on pressure drop is a property of Newtonian fluids. Again, the Newtonian assumption is only applicable to larger diameter blood vessels. The identical results obtained for branches “2” and “3” validate the complicated programming logic and units conversion code.

Figure 3 Calculated results appear in gray boxes (symmetrical bifurcation data are used to validate computed results for required symmetric solutions).

2.2.2 A Second Calculation

Now we consider for purposes of validation the calculation in Figure 4. Here, the branches in “2” and “3” are again identical, but longer and narrower than that in “1.” The pressure drop in “1” is now 2.12684 – 2.11938 = 0.00746 psi, while that in “2” and “3” is 2.11938 – 2 = 0.11938 psi. Even though the flow rate in “1” is twice that of “2” or “3,” the pressure drop is much smaller, because the radius in “1” is twice as large. The total pressure drop from inlet to outlet is 2.12684 – 2 = 0.12684 psi. Another Newtonian flow property is worth mentioning. If the viscosity is halved to 5 cp, the total pressure drop becomes 2.06342 – 2 = 0.06342 psi, which is half of 0.12684 psi (software screen not shown). In other words, a decrease in viscosity leads to a proportional decrease in pressure drop. The “scalability” properties demonstrated above do not apply to non-Newtonian fluids. While different rheology models may also be required as suggested by empirical data, in general for such fluids, computational methods are required when geometric rate or pressure data is altered.

Figure 4 Branch “1” different from “2” and “3,” which are identical.

2.3 Theory – Two Uneven Bifurcated Arteries with P_i Specified

The above applies when Q₁ is the main input. We showed earlier how if P_i, P_o,2 and P_o,3 are known, the junction pressure follows from P_a = (R₁⁴ P_i /L₁ + R₂⁴ P_o,2 /L₂ + R₃⁴ P_o,3 /L₃)/(R₁⁴ /L₁ + R₂⁴ /L₂ + R₃⁴ /L₃). Then, the inlet flow rate follows from $\text{Q}_{1}=\unicode[Times]{x03C0} {\text{R}_{2}}^{4}\left(\text{P}_ \mathcal a-\text{P}_{o, 2}\right) /\left(8 \unicode[Times]{x03BC} \text{L}_{2}\right)+\unicode[Times]{x03C0} \text{R}_{3}{ }^{4}\left(\text{P}_ \mathcal a-\text{P}_{o, 3}\right) /\left(8 \unicode[Times]{x03BC} \text{L}_{3}\right)$. With these quantities available, an analysis similar to the above is straightforward.

2.4 Software – CODE-2 for Two Uneven Bifurcated Arteries with P_i Specified

We now discuss the software implementation where P_i is given, referred to as “CODE-2.” Obviously, CODE-1 and CODE-2 are complementary. They were developed separately to keep the programming logic simple. Figure 5 shows both side-by-side while they are creating consistent numerical outputs to the fourth decimal place.

Figure 5 Complementary CODE-1 and CODE-2 implementations.

2.5 Theory – Complicated Network Flows with Chained Bifurcations

Here we consider a more complicated network with additional bifurcations, as shown in Figure 6. It is obviously not the most general, nor was it intended to be. Our objective is to demonstrate the algebraic manipulations needed to handle networks of greater complexity, so readers can develop custom models as needed. For clarity, we will drop all “outlet” subscripts, e.g., P_o,2 becomes P₂.

Figure 6 More complicated flow network.

As before, we start with the Hagen-Poiseuille law, as given in Equations 2.13 and 2.14. The parameter “A” is introduced for simplicity, to eliminate an error-prone “ π/(8μ) .” Whereas we invoked “Q₁ = Q₂ + Q₃” before, we now perform the analogous operation twice, as in Equations 2.15 and 2.16 – that is, once for each network, focusing on the junctions a and b.

\[ \text{Q}=\unicode[Times]{x03C0} \frac{\text{R}^{4}\left(\text{P}_{i}-\text{P}_{o}\right)}{8 \unicode[Times]{x03BC} L} \tag{2.13} \]

\[ \text{Q}=\frac{\text{AR}^{4}\left(\text{P}_{i}-\text{P}_{o}\right)}{\text{L}} \text { where } \text{A}=\frac{\unicode[Times]{x03C0}}{8 \unicode[Times]{x03BC}} \tag{2.14} \]

\[ \mathrm{Q}_{1 a}=\mathrm{Q}_{a 2}+\mathrm{Q}_{a 3}+\mathrm{Q}_{a b} \tag{2.15} \]

\[ \mathrm{Q}_{a b}=\mathrm{Q}_{b 4}+\mathrm{Q}_{b 5}+\mathrm{Q}_{b 6}+\mathrm{Q}_{b 7} \tag{2.16} \]

We next substitute Equation 2.14 into the right side of Equation 2.15, and Equation 2.14 into all terms in Equation 2.16. This leads to

\[ \mathrm{Q}_{1 a}=\frac{\mathrm{AR}_{a 2}{ }^{4}\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{2}\right)}{\mathrm{L}_{a 2}}+\frac{\mathrm{AR}_{a 3}{ }^{4}\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{3}\right)}{\mathrm{L}_{a 3}}+\mathrm{AR}_{a b}{ }^{4} \frac{\mathrm{P}_ \mathcal a-\mathrm{P}_{b}}{\mathrm{L}_{a b}} \tag{2.17} \]

\[ \frac{\mathrm{R}_{a b}{ }^{4}\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{b}\right)}{\mathrm{L}_{a b}}=\frac{\mathrm{R}_{b 4}{ }^{4}\left(\mathrm{P}_{b}-\mathrm{P}_{4}\right)}{\mathrm{L}_{b 4}}+\frac{\mathrm{R}_{b 5}{ }^{4}\left(\mathrm{P}_{b}-\mathrm{P}_{5}\right)}{\mathrm{L}_{b 5}}+\frac{\mathrm{R}_{b 6}{ }^{4}\left(\mathrm{P}_{b}-\mathrm{P}_{6}\right)}{\mathrm{L}_{b 6}}+\frac{\mathrm{R}_{b 7}{ }^{4}\left(\mathrm{P}_{b}-\mathrm{P}_{7}\right)}{\mathrm{L}_{b 7}} \tag{2.18} \]

Interestingly, the above can be rewritten as “two equations for two unknowns” in the elementary algebra sense. The unknowns are the junction pressures P_a and P_b, as shown in Equations 2.19 and 2.20, with coefficients B, C, D, E, F and G as in Equations 2.24 – 2.26. The solutions are given in Equations 2.27 and 2.28.

\[ \mathbf{B}\text{P}_ \mathcal a+\mathbf{C }\text{P}_{b}=\mathbf{D} \tag{2.19} \]

\[ \mathbf{E}\text{P}_ \mathcal a+\mathbf{F}\text{P}_{b}=\mathbf{G} \tag{2.20} \]

\[ \mathbf{B}=\mathrm{A}\left( \frac{\mathrm{{R}_{a 2}}^{4}}{\mathrm{L}_{a 2}}+ \frac{\mathrm{{R}_{a 3}}^{4}}{\mathrm{L}_{a 3}}+ \frac{\mathrm{{R}_{a \textit{b}}}^{4}}{\mathrm{L}_{a \textit{b}}}\right) \tag{2.21} \]

\[ \mathbf{C}=-\mathrm{A} \frac{\mathrm{{R}_{a \textit{b}}}^{4}}{\mathrm{L}_{a \textit{b}}} \tag{2.22} \]

\[ \mathbf{D}=\mathrm{Q}_{1 a}+\frac{\mathrm{{AR}_{a 2}}^{4} \mathrm{P}_{2}}{\mathrm{L}_{a 2}}+\frac{\mathrm{{AR}_{a 3}}^{4} \mathrm{P}_{3}}{\mathrm{L}_{a 3}} \tag{2.23} \]

\[ \mathbf{E}=\frac{\mathrm{{R}_{a \textit{b}}}^{4}}{\mathrm{L}_{a \textit{b}}} \tag{2.24} \]

\[ \mathbf{F}=-\left(\frac{\mathrm{{R}_{a \textit{b}}}^{4}}{\mathrm{L}_{a \textit{b}}}+ \frac{{\mathrm{R}_{b 4}}^{4}}{\mathrm{L}_{b 4}}+ \frac{{\mathrm{R}_{b 5}}^{4}}{\mathrm{L}_{b 5}}+ \frac{{\mathrm{R}_{b 6}}^{4}}{\mathrm{L}_{b 6}}+ \frac{{\mathrm{R}_{b 7}}^{4}}{\mathrm{L}_{b 7}}\right) \tag{2.25} \]

\[ \mathbf{G}=-\left(\frac{\mathrm{R}_{b 4}{ }^{4} \mathrm{P}_{4}}{\mathrm{L}_{b 4}}+\frac{\mathrm{R}_{b 5}{ }^{4} \mathrm{P}_{5}}{\mathrm{L}_{b 5}}+\frac{\mathrm{R}_{b 6}{ }^{4} \mathrm{P}_{6}}{\mathrm{L}_{b 6}}+\frac{\mathrm{R}_{b 7}{ }^{4} \mathrm{P}_{7}}{\mathrm{L}_{b 7}}\right) \tag{2.26} \]

\[ \mathrm{P}_ \mathcal a=\frac{(\mathbf{D F}-\mathbf{C G})}{\mathbf{B F}-\mathbf{C E}} \tag{2.27} \]

\[ \mathrm{P}_{b}=\frac{(\mathbf{B G}-\mathbf{D E})}{\mathbf{B F}-\mathbf{C E}} \tag{2.28} \]

The final steps involve post-processing for inlet pressure, volume flow rates and shear stresses. First, from Equation 2.14, the relationship “Q = AR⁴(P_i – P_o)/L” suggests that Q₁ = AR₁⁴(P_i – P_a)/L₁. Since Q₁ is prescribed and P_a is known from above, this yields

\[ \text{P}_{i}=\text{P}_ \mathcal a+\frac{\text{Q}_{1} \text{L}_{1}}{\text{A} \text{R}_{1}{ }^{4}} \tag{2.29} \]

It suffices to consider a typical branch from each of junctions “a” and “b.” For example, from “a,” we have

\[ \mathrm{Q}_{2}=\mathrm{AR}_{2}{ }^{4} \frac{\mathrm{P}_ \mathcal a-\mathrm{P}_{2}}{\mathrm{L}_{2}} \tag{2.30} \]

\[ \unicode[Times]{x03C4}_2=\frac{\left(\text{P}_ \mathcal a-\text{P}_{2}\right) \text{R}_{2}}{2 \text{L}_{2}} \tag{2.31} \]

while from “b,” we have

\[ \mathrm{Q}_{6}=\mathrm{{AR}_{6}}^{4} \frac{\mathrm{P}_{b}-\mathrm{P}_{6}}{\mathrm{L}_{6}} \tag{2.32} \]

\[ \unicode[Times]{x03C4}_{6}=\frac{\left(\text{P}_{b}-\text{P}_{6}\right) \text{R}_{6}}{2 \text{L}_{6}} \tag{2.33} \]

Software is straightforwardly developed from the above for Figure 6 and similar networks. We leave this to the interested reader.

2.6 Network with Arbitrary Number of Bifurcations

Here we consider an earlier problem where Q₁ is given, and the junction pressure P_a and inlet pressure P_i are required, but now an arbitrary number of bifurcations (having different radii and lengths) with n = 2, 3, 4 . . . N is allowed. Volume flow rates and shear stresses in all sub-branches will be derived. More than likely, the most popular application will be the earlier two-branch model whose solution was worked out in detail. However, the present generalization will be useful for predictions in more complicated networks.

We begin by expressing mass conservation Q₁ = Q₂ + Q₃ + . . . + Q_N in summation notation, where the sum ∑ is taken over n = 2, 3, 4 . . . N, that is,

\[ \mathrm{Q}_{1}=\left\{\frac{\unicode[Times]{x03C0}}{8 \unicode[Times]{x03BC}}\right\} \sum \mathrm{{R}_{\mathrm{n}}}^{4} \frac{\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{\mathrm{n}}\right)}{\mathrm{L}_{\mathrm{n}}} \tag{2.34} \]

The constant P_a can be factored out and solved as

\[ \text{P}_ \mathcal a=\frac{\left(8 \unicode[Times]{x03BC} \frac{\text{Q}_{1}}{\unicode[Times]{x03C0}}+\sum \frac{\text{R}_{\text{n}}{ }^{4} \text{P}_{\text{n}}}{\text{L}_{\text{n}}}\right)}{\sum\left(\frac{\text{R}_{\text{n}}{ }^{4}}{\text{L}_{\text{n}}}\right)} \tag{2.35} \]

If we now re-express Equation 2.34 in the expanded form R₁⁴(P_i – P_a)/L₁ = ∑ R_n⁴ (P_a – P_n) /L_n, we find the inlet pressure as

\[ \mathrm{P}_{\mathrm{i}}=\mathrm{P}_ \mathcal a\left[1+\left(\frac{\mathrm{L}_{1}}{\mathrm{R}_{1}{ }^{4}}\right) \sum \frac{\mathrm{R}_{\mathrm{n}}{ }^{4}}{\mathrm{L}_{\mathrm{n}}}\right]-\left(\frac{\mathrm{L}_{1}}{\mathrm{R}_{1}{ }^{4}}\right) \sum \frac{\mathrm{R}_{\mathrm{n}}{ }^{4} \mathrm{P}_{\mathrm{n}}}{\mathrm{L}_{\mathrm{n}}} \tag{2.36} \]

where the P_a on the right side is known from Equation 2.35. Then, the flow rates and viscous shear stresses below are obtained.

\[ \mathrm{Q}_{\mathrm{n}}=\unicode[Times]{x03C0} \mathrm{R}_{\mathrm{n}}{ }^{4} \frac{\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{\mathrm{n}}\right)}{8 \unicode[Times]{x03BC} \mathrm{L}_{\mathrm{n}}} \tag{2.37} \]

\[ \unicode[Times]{x03C4}_1=\frac{\left(\text{P}_{i}-\text{P}_ \mathcal a\right) \text{R}_{1}}{2 \text{L}_{1}} \tag{2.38} \]

\[ \unicode[Times]{x03C4}_{\mathrm{n}}=\frac{\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{\mathrm{n}}\right) \mathrm{R}_{\mathrm{n}}}{2 \mathrm{L}_{\mathrm{n}}} \tag{2.39} \]

Extensions to networks more complicated than that of Figure 6 are straightforward. The “summation strategy” above would apply to a sequence “a, b, c, d . . .” of junction point problems, and coupled linear algebraic equations for P_a, P_b, P_c, P_d and so on, would be solved. Then, inlet pressure is derived, followed by results for all branch volume flow rates and viscous shear stresses.

2.7 Bifurcated Newtonian Flow in Noncircular Clogged Blood Vessels

We have so far considered bifurcation models where the flow cross-sections are perfectly circular. This applies to clean circular pipes and blood vessels. However, these may clog, resulting in arbitrarily formed duct geometries. Typical examples are given in Figure 7a and Figure 7b. The general clog analysis model of Chin [12] first maps irregular domains into simple rectangular spaces associated with boundary-conforming curvilinear grids. The transformed equations are then solved in this domain. Figure 7c displays time-lapse results obtained for a blood clotting simulation where plaque buildup is limited by transient viscous shear stress and dynamic erosive actions (the pseudo-transient analysis ignores fluid density). The reader is referred to the cited reference, which deals primarily with single conduit flows, for details. For the purposes of this paper, it suffices to note that in the case of clogged Newtonian flows, the classical Hagen-Poiseuille formula is modified by a multiplicative correction factor and an effective radius can be defined. Then, the results obtained above immediately apply.

Figure 7 a) Irregular clog distributions and debris buildup in pipelines and conduits. b) Bifurcated artery with multiple clogged flow branches. c) Time-lapse, pseudo-transient clog formation in non-Newtonian fluid.

It is important to understand that relationships like “flow rate versus pressure drop” and shear stress details depend more on the details of the fluid and clog interface, and to a much lesser extent, on “clogged versus unclogged area ratios.” The required clog modifications are simple. Suppose, for given numerical values dp/dz* and μ∗ we fix the clog geometry and calculate the corresponding flow rate Q* using the gridding method of Chin [13]. As shown in that reference, Q* = - GF* (dp/dz*/μ*) applies to Newtonian flows where GF* is a geometric factor representative of the cross-section. We calculate the value GF* = - Q*/(dp/dz*/μ*). For a Hagen-Poiseuille flow, we know that the flow rate Q = - {πR⁴ /8}{dp/dz/μ }. Thus, we can introduce an effective radius R_eff = (8 GF*/π )^1/4 (assuming like values of "dp/dz/μ ") for each clogged flow element in the analytical bifurcation results above. Note that shear stress formulas like Equations 2.6-2.8, which assume axisymmetry, do not apply. However, as seen from Figure 7c, non-uniform viscous stresses are available everywhere in the flow domain.

3. Non-Newtonian Flow in Circular Conduits and Networks

In this section, we consider non-Newtonian flows in circular conduits and networks formed from such sections. Because nonlinear relations connect flow rates to pressure differentials, simple analytical solutions are not possible and recourse is made to numerical methods. We first illustrate our approach for the simplest “Power Law” model. Then, corresponding models are developed for Bingham Plastic and Herschel-Bulkley and detailed explanations are offered.

3.1 Power Law Fluids

In the simplest case, the viscous shear stress $\unicode[Times]{x03C4}$ and the axial velocity profile U(r) are related by $\unicode[Times]{x03C4}$(r) = k {dU/dr}ⁿ where r is the radial coordinate, “n” is a dimensionless exponent and "k" is the “consistency index.” This leads to the volume flow rate

\[ \text{Q}=\left\{\text{n} \unicode[Times]{x03C0} \frac{\text{R}^{3}}{3 \text{n+1}}\right\}\left[\frac{\text{R}\left(\text{P}_{i}-\text{P}_{o}\right)}{2 \text{k}\text{L} }\right]^{1 / \text{n}} \tag{3.1} \]

The Newtonian Hagen-Poiseuille formula is obtained by setting n = 1 and the consistency factor to k =μ , where μ is the constant viscosity. This leads to the $\text{Q}_{}=\unicode[Times]{x03C0} \text{R}^{4}\left(\text{P}_{i}-\text{P}_{o}\right) /\left(8 \unicode[Times]{x03BC} \text{L}_{}\right)$. The assumptions for the Power Law model are similar. They are: (a) a rigid circular cross-section of radius R, (b) a straight conduit with length L satisfying L >> R, (c) constant values of k and n, (d) end conditions with P_i and P_o respectively being inlet and outlet pressures, and (e) steady-state flow. Our convention assumes that the volume flow rate Q > 0 if P_i > P_o. Note that fluid density does not appear because steady-state, laminar (non-turbulent) conditions are assumed – density, however, will appear in transient and compressible flow applications.

We consider the two-branch bifurcation shown in Figure 8. Conservation of mass requires that Q₁ = Q₂ + Q₃ where Q₁ is specified so that we may write

\[ \mathrm{Q}_{1}=\left\{\mathrm{n} \unicode[Times]{x03C0} \frac{\mathrm{R}_{2}{ }^{3}}{3 \mathrm{n}+1}\right\}\left[\frac{\mathrm{R}_{2}\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{\mathrm{o}, 2}\right)}{2 \mathrm{~kL}_{2}}\right]^{1 / \mathrm{n}}+\left\{\mathrm{n} \unicode[Times]{x03C0} \frac{\mathrm{R}_{3}{ }^{3}}{3 \mathrm{n}+1}\right\}\left[\frac{\mathrm{R}_{3}\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{\mathrm{o}, 3}\right)}{2 \mathrm{~kL}_{3}}\right]^{1 / \mathrm{n}} \tag{3.2} \]

Figure 8 Simple two-branch bifurcation.

3.1.1 Iterative “Half-step” Solution for P_a and P_i

Earlier, we solved for the junction pressure P_a initially, and using it, derived expressions for Q₂ and Q₃. The algebra was simple and lengthy manipulations were performed. Here, the exponent “n” is general, possibly taking on integer and fractional values, so that an expression for P_a could not be obtained. In fact, Equation 3.2 is nonlinear. One might attempt to use Newton-Raphson schemes popular in solving high-order polynomial equations, but these are slowly convergent and rely on the availability of an initially close solution. Any method we devise must be robust and yield solutions with minimal user interaction while furthermore applying to other rheological models used for our problems.

When Q₁, P_o,2 and P_o,3 are given, we would expect Q₁ to increase as P_a increases. This monotonic behavior can be put to use in deriving the needed solution. We simply start with a guess, the best being the greater of P_o,2 and P_o,3, to ensure that the quantities in the square brackets of Equation 3.2 are both positive. Then, as P_a is successively increased by given ∆P increments, the right hand side “RHS” of Equation 3.2 will increase until the ratio “RHS/Q_1,prescribed” exceeds unity. This means that we have over-predicted both P_a and Q₁. To compensate for this, we return to our “previous P_a” and guess using a reduced ΔP/2 instead. If we still over-predict, the iterative process is repeated, with this recursive “half step” operation continued until a user-specified error bound is satisfied.

Once we have obtained P_a, we turn to Equation 3.1 to infer the form for volume flow rate expressions in “1” and the two downstream branches “2” and “3,” that is

\[ \text{Q}_{1}=\left\{\text{n} \unicode[Times]{x03C0} \frac{\text{R}_{1}{ }^{3}}{3 \text{n+1}}\right\}\left[\text{R}_{1} \frac{\left(\text{P}_{i}-\text{P}_ \mathcal a\right)}{2 \text{k} \text{L}_{1}}\right]^{1 / \text{n}} \tag{3.3} \]

\[ \text{Q}_{2}=\left\{\text{n} \unicode[Times]{x03C0} \frac{\text{R}_{2}{ }^{3}}{3 \text{n+1}}\right\}\left[\text{R}_{2} \frac{\text{P}_ \mathcal a-\text{P}_{o,2}}{2 \text{k} \text{L}_{2}}\right]^{1 / \text{n}} \tag{3.4} \]

\[ \text{Q}_{3}=\left\{\text{n} \unicode[Times]{x03C0} \frac{\text{R}_{3}{ }^{3}}{3 \text{n+1}}\right\}\left[\text{R}_{3} \frac{\text{P}_ \mathcal a-\text{P}_{o,3}}{2 \text{k} \text{L}_{3}}\right]^{1 / \text{n}} \tag{3.5} \]

Equation 3.3 can be solved to produce the inlet pressure P_i for branch “1,” with the result that

\[ \text{P}_{i}=\text{P}_ \mathcal a+\left(\frac{2 \text{k} \text{L}_{1}}{\text{R}_{1}}\right)\left[\frac{(3 \text{n}+1) \text{Q}_{1}}{\text{n} \unicode[Times]{x03C0} \text{R}_{1}{ }^{3}}\right]^{n} \tag{3.6} \]

3.1.2 Shear Stress

For Power Law pipe flows, well known results for laminar viscous shear stress can be used. In our convention used, $\unicode[Times]{x03C4}$= k (– dU/dr) ⁿ where the velocity profile satisfies the equation U(r) = ((P_i – P_a)/(2kL₁)) ^1/n [n/(n+1)] (R₁ ^(n+1)/n – r ^(n+1)/n). This leads to the derivative relation (– dU/dr)ⁿ = ((P_i – P_a)/(2kL₁)) r, which in turn yields the formula $\unicode[Times]{x03C4}$= k(– dU/dr)ⁿ = (P_i – P_a) R₁/(2L₁) > 0 at r = R₁. Similar results are obtained at R₂ and R₃, so that,

\[ \unicode[Times]{x03C4}_1=\left(\mathrm{P}_{\mathrm{i}}-\mathrm{P}_ \mathcal a\right) \frac{\mathrm{R}_{1}}{2 \mathrm{L}_{1}} \tag{3.7} \]

\[ \unicode[Times]{x03C4}_2=\left(\mathrm{P}_{\mathrm {\mathcal a}}-\mathrm{P}_{o,2}\right) \frac{\mathrm{R}_{2}}{2 \mathrm{L}_{2}} \tag{3.8} \]

\[ \unicode[Times]{x03C4}_3=\left(\mathrm{P}_{\mathrm {\mathcal a}}-\mathrm{P}_{o,3}\right) \frac{\mathrm{R}_{3}}{2 \mathrm{L}_{3}} \tag{3.9} \]

These results are identical to those for Newtonian fluids. In fact, they follow elegantly from global considerations balancing pressure force and wall shear effects. However, for complicated geometries containing general clogs, the details of the clog and fluid interface enter and the above formulas do not apply. These problems are considered in Chin [12] and Chin and Chin [13].

3.1.3 Typical Parameters

The solution process for P_a outlined earlier can be written in Fortran, as shown in Figure 9, which can be translated into other computer languages or spreadsheet formats. Five “canned” examples are provided, which are described in the following pages. We note that the basic physical units LBF, IN and SEC for pound force, inch and second are used in our analyses. Pressure, for example, is then represented by PSI or LBF/IN². For water at room temperature and pressure, the Newtonian rheology parameters are n = 1 and about k =μ = 0.0000001465 psi sec. This flow is always used in our basic validations. Examples related to a blood flow application are described in Chin and Chin [13]. We will follow the assumptions in that example – an average consistency index of k = 17.0 mPa sec ⁿ that has been determined by medical investigators for healthy control subjects. In our units, this is 17 (0.0000001450 psi sec ⁿ) or 0.000002465 psi sec ⁿ – a bit less than 17 cp, if n were equal to unity. They also found n values slightly greater than 0.7, supporting their observation that blood mixtures are Power Law in nature. Others have assumed viscosities in the 10-15 cp range, with similar findings for n. To assess the robustness of our algorithms, we considered high and low sides of n and k values. Note that Power Law models neglect yield stresses, which may be important in flow blockage and clotting. When these are important, the Power Law model becomes a Herschel-Bulkley fluid, which will be considered later. We next turn to the algorithm BIFURC-4 for junction pressure P_a, shown in its entirety in Figure 9.

Figure 9 Power Law solution for junction pressure.

In the above code, the Power Law solution for junction pressure P_a (psi) is provided where Q₁ (cc/s) is given, with iterations defined by simple “DO 100” loop. We now describe the iterative “half-step” solutions obtained for P_a. In Newtonian flow Example 1 immediately below, the assumed computational parameters are identical to those used in the exact solution underlying Figure 3. The final solution, highlighted in bold font, yields = 2.00373 psi, which also appears in Figure 3 for our assumed Q₁ = 100 cc/s (note how the “test Q” almost matches the given 100 cc/s). In the examples below, “N” is a Fortran “counter” which records the “guess number” for P_a, whose value appears immediately to its right. As N increases, the ratio monotonically increases, with the calculations terminating once the RATIO reaches unity to within a user-prescribed error bound. The computational histories are robust and stable. Despite the large N values used for Examples 4 and 5, the desk time involved for a Windows i5 computer is much less than one second.

C:\bifurc-4

Example 1-5?   Enter:   1

   N           PA                  QTEST           QCCPS           RATIO

369   0.200370E+01   0.9905E+02   0.1000E+03   0.9904621

370   0.200371E+01   0.9931E+02   0.1000E+03   0.9931462

371   0.200372E+01   0.9958E+02   0.1000E+03   0.9958304

372   0.200373E+01   0.9985E+02   0.1000E+03   0.9985145

Stop - Program terminated.

Example 1-5?   Enter:   2

   N           PA                  QTEST           QCCPS           RATIO

.

1146   0.201148E+01   0.9942E+02   0.1000E+03   0.9941543

1147   0.201149E+01   0.9954E+02   0.1000E+03   0.9953938

1148   0.201150E+01   0.9966E+02   0.1000E+03   0.9966338

1149   0.201151E+01   0.9979E+02   0.1000E+03   0.9978743

1150   0.201152E+01   0.9991E+02   0.1000E+03   0.9991152

Stop - Program terminated.

Example 1-5?   Enter:   3

   N           PA                  QTEST           QCCPS           RATIO

.

12008   0.201145E+01   0.9912E+02   0.1000E+03   0.9912050

12009   0.201145E+01   0.9913E+02   0.1000E+03   0.9913229

12010   0.201145E+01   0.9914E+02   0.1000E+03   0.9914407

.

.

12060   0.201150E+01   0.9973E+02   0.1000E+03   0.9973426

12070   0.201151E+01   0.9985E+02   0.1000E+03   0.9985242

12080   0.201152E+01   0.9997E+02   0.1000E+03   0.9997062

12081   0.201152E+01   0.9998E+02   0.1000E+03   0.9998245

12082   0.201152E+01   0.9999E+02   0.1000E+03   0.9999428

Stop - Program terminated.

Example 1-5?   Enter:   4

   N           PA                  QTEST           QCCPS           RATIO

19227   0.219253E+01   0.9901E+02   0.1000E+03   0.9900671

19228   0.219254E+01   0.9901E+02   0.1000E+03   0.9901407

19229   0.219255E+01   0.9902E+02   0.1000E+03   0.9902143

19230   0.219256E+01   0.9903E+02   0.1000E+03   0.9902877

.

19330   0.219356E+01   0.9977E+02   0.1000E+03   0.9976528

19340   0.219366E+01   0.9984E+02   0.1000E+03   0.9983900

19350   0.219376E+01   0.9991E+02   0.1000E+03   0.9991276

19360   0.219386E+01   0.9999E+02   0.1000E+03   0.9998654

19361   0.219387E+01   0.9999E+02   0.1000E+03   0.9999391

Stop - Program terminated.

Example 1-5?   Enter:   5

   N           PA                  QTEST           QCCPS           RATIO

.

8780   0.208792E+01   0.9991E+02   0.1000E+03   0.9991337

8781   0.208793E+01   0.9994E+02   0.1000E+03   0.9993613

8782   0.208794E+01   0.9996E+02   0.1000E+03   0.9995889

8783   0.208795E+01   0.9998E+02   0.1000E+03   0.9998167

Stop - Program terminated.

Having successfully determined the junction pressure P_a, we now evaluate Equations 3.3 – 3.9 for inlet pressure, branch volume flow rates and wall viscous shear stresses. In our BIFURC-5 code (not shown), we added the Fortran logic below and re-compiled the source code.

105  PPSI1 = PA + (2.*K*L1/R1)*

1      ((3.*EN+1.)*QIN3PS/(EN*PI*R1**3.))**EN

Q1 = QCCPS

Q2 = 16.38706*(EN*PI*R2**3./(3.*EN+1.))*

1      (R2*(PA-P2)/(2.*K*L2))**EXPON

Q3 = 16.38706*(EN*PI*R3**3./(3.*EN+1.))*

1      (R3*(PA-P3)/(2.*K*L3))**EXPON

TAU1 = (PPSI1 - PA)*R1/(2.*L1)

TAU2 = (PA - P2)*R2/(2.*L2)

TAU3 = (PA - P3)*R3/(2.*L3)

We now re-run the prior Examples 1-5. In Example 1 immediately below, we had considered a purely Newtonian flow, with a 10 cp viscosity, previously given in Figure 3.

C:\>BIFURC-5

Example 1-5? Enter: 1

   N           PA                  QTEST           QCCPS           RATIO

370   0.200371E+01   0.9931E+02   0.1000E+03   0.9931462

371   0.200372E+01   0.9958E+02   0.1000E+03   0.9958304

372   0.200373E+01   0.9985E+02   0.1000E+03   0.9985145

P1 = 2.0111911 (PSI)

Q1 = 100.00, Q2 = 49.99, Q3 = 49.99 (CC/S)

TAU1 = 0.00018653,TAU2 = 0.00009325, TAU3 = 0.00009325 (PSI)

Stop - Program terminated.

The results are excellent. Earlier results were based on exact, closed form, analytical solutions using the Newtonian model, while here, they follow from a “half-step” iterative numerical scheme for Power Law fluids. The exact inlet pressure was earlier seen to be 2.01119 psi, and the iterated value obtained here is also 2.01119 psi. Also, the main flow is correctly divided into two branches having equal 49.99 cc/s (as compared to 50 cc/s) rates. Finally, the equal shear stresses in branches “2” and “3” agree with exact solutions to the fourth decimal place. Computed results for Examples 2-5 are given below.

Example 1-5? Enter: 2

   N           PA                  QTEST           QCCPS           RATIO

.

.

1148   0.201150E+01   0.9966E+02   0.1000E+03   0.9966338

1149   0.201151E+01   0.9979E+02   0.1000E+03   0.9978743

1150   0.201152E+01   0.9991E+02   0.1000E+03   0.9991152

P1 = 2.0302393 (PSI)

Q1 = 100.00, Q2 = 49.99, Q3 = 49.99 (CC/S)

TAU1 = 0.00046797, TAU2 = 0.00028802, TAU3 = 0.00028802(PSI)

Stop - Program terminated.

Example 1-5? Enter: 3

   N           PA                  QTEST           QCCPS           RATIO

.

12079   0.201152E+01   0.9996E+02   0.1000E+03   0.9995880

12080   0.201152E+01   0.9997E+02   0.1000E+03   0.9997062

12081   0.201152E+01   0.9998E+02   0.1000E+03   0.9998245

12082   0.201152E+01   0.9999E+02   0.1000E+03   0.9999428

P1 = 2.0302415 (PSI)

Q1 = 100.00, Q2 = 50.00, Q3 = 50.00 (CC/S)

TAU1 = 0.00046797, TAU2 = 0.00028807, TAU3 = 0.00028807(PSI)

Stop - Program terminated.

Example 1-5? Enter: 4

   N           PA                  QTEST           QCCPS           RATIO

.

19358   0.219384E+01   0.9997E+02   0.1000E+03   0.9997178

19359   0.219385E+01   0.9998E+02   0.1000E+03   0.9997916

19360   0.219386E+01   0.9999E+02   0.1000E+03   0.9998654

19361   0.219387E+01   0.9999E+02   0.1000E+03   0.9999391

P1 = 2.5088389 (PSI)

Q1 = 100.00, Q2 = 50.00, Q3 = 50.00 (CC/S)

TAU1 = 0.00787402, TAU2 = 0.00484695, TAU3 = 0.00484695(PSI)

Stop - Program terminated.

Example 1-5? Enter: 5

   N           PA                  QTEST           QCCPS           RATIO

.

.

8779   0.208791E+01   0.9989E+02   0.1000E+03   0.9989061

8780   0.208792E+01   0.9991E+02   0.1000E+03   0.9991337

8781   0.208793E+01   0.9994E+02   0.1000E+03   0.9993613

8782   0.208794E+01   0.9996E+02   0.1000E+03   0.9995889

8783   0.208795E+01   0.9998E+02   0.1000E+03   0.9998167

P1 = 2.2123446 (PSI)

Q1 = 100.00, Q2 = 50.00, Q3 = 50.00 (CC/S)

TAU1 = 0.00310976, TAU2 = 0.00219886, TAU3 = 0.00219886(PSI)

Stop - Program terminated.

Although there are no Power Law solutions with which to benchmark the above results, the calculated inlet pressures and shear stresses appear to be correct in order-of-magnitude as suggested by our Newtonian results. We note that the iterations are rapidly convergent and require less than a second on Windows computers.

3.2 Herschel-Bulkley Fluids

In developing math models and numerical solutions, it is important to “start simple” and work towards more complicated problems incrementally. Validation is key at every step of the process. For example, we had started with simple Newtonian flow, first considering a two-branch bifurcation with Q₁ specified initially, and then strategies for additional branches and Power Law fluids were considered.

3.2.1 Analytical and Numerical Approach

We had first dealt with nonlinear Power Law fluids. Because exact, closed form, analytical solutions similar to those found in Newtonian flow are not possible and Newton-Raphson root finder methods are not ideal, we devised a straightforward “half-step” solution that takes advantage of the monotonic relation between pressure gradient and flow rate. This method was validated against exact analytical results using a Newtonian example – such fluids, again, represent one limit of the more general Power Law model. We now consider Herschel-Bulkley fluids, which extend Power Law descriptions to include yield stress. The mathematics are more involved, so we will similarly validate new results against earlier Newtonian and Power Law calculations.

Yield Stress Modeling. The classical Herschel-Bulkley pipe flow model combines Power Law with yield stress ($\unicode[Times]{x03C4}$_y) characteristics, with the result that

\[ \unicode[Times]{x03C4}=\unicode[Times]{x03C4}_{y}+ \text{k}\left(-\frac{\text{dU}}{\text{dr}}\right)^{n} \tag{3.10} \]

\[ \mathrm{U}(\mathrm{r})=\mathrm{k}^{-1 / \mathrm{n}}\left(\frac{\Delta \mathrm{P}}{2 \mathrm{L}}\right)^{-1}\left\{\frac{\mathrm{n}}{\mathrm{n}+1}\right\} \]

\[ \times\left[\left(\mathrm{R} \frac{\Delta \mathrm{P}}{2 \mathrm{L}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right)^{(\mathrm{n}+1) / \mathrm{n}}-\left(\mathrm{r} \frac{\Delta \mathrm{P}}{2 \mathrm{L}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right)^{(\mathrm{n}+1) / \mathrm{n}}\right], \mathrm{R}_{\mathrm{p}} \leq \mathrm{r} \leq \mathrm{R} \tag{3.11} \]

\[ \text{U}(\mathrm{r})=\text{k}^{-1 / \text{n}}\left(\frac{ \Delta \text{P}}{2\text{L}}\right)^{-1}\left\{\frac{\text{n}}{\text{n}+1}\right\} \]

\[ \times\left[\left(\mathrm{R} \frac{\Delta \mathrm{P}}{2 \mathrm{L}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right)^{(\mathrm{n}+1) / \mathrm{n}}-\left(\mathrm{R}_{\mathrm{p}} \frac{\Delta \mathrm{P}}{2 \mathrm{L}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right)^{(\mathrm{n}+1) / \mathrm{n}}\right], 0 \leq \mathrm{r} \leq \mathrm{R}_{\mathrm{p}} \tag{3.12} \]

\[ \text{R}_{\text{p}}=2 \unicode[Times]{x03C4}_{y} \frac{\text{L}}{ \Delta \text{P}} \tag{3.13} \]

\[ \frac{\text{Q}}{\unicode[Times]{x03C0} \text{R}^{3}}=\text{k}^{-1 / \text{n}}\left(\text{R} \frac{ \Delta \text{P}}{2\text{L}}\right)^{-3}\left(\text{R} \frac{ \Delta \text{P}}{2\text{L}}-\unicode[Times]{x03C4}_{y}\right)^{ \text{(n+1)} / \text{n}} \]

\[ \times\left[\left(\text{R} \frac{ \Delta \text{P}}{2\text{L}}-\unicode[Times]{x03C4}_{y}\right)^{2} \frac{\text{n}}{3 \text{n}+1}+2 \unicode[Times]{x03C4}_{y}\left(\text{R} \frac{ \Delta \text{P}}{2\text{L}}-\unicode[Times]{x03C4}_{y}\right) \frac{\text{n}}{2 \text{n}+1}+\unicode[Times]{x03C4}_{y}{ }^{2} \frac{\text{n}}{\text{n}+1}\right] \tag{3.14} \]

Here, the “plug radius” R_p defines the outer extent of the cylinder that moves as a solid body, residing at the center of the pipe. Because flows with less than and greater than this radius are different, two solutions for the axial velocity U(r) appear above. Our usual starting point, namely, mass conservation assuming “Q₁ = Q₂ + Q₃,” requires first a volume integration of U(r) to produce Q(r), and second, the explicit representation

\[ \mathrm{Q}_{1}=\unicode[Times]{x03C0} \mathrm{R}_{2}{ }^{3} \mathrm{k}^{-1 / \mathrm{n}}\left(\frac{\mathrm{R}_{2}\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{2}\right)}{2 \mathrm{L}_{2}}\right)^{-3}\left(\frac{\mathrm{R}_{2}\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{2}\right)}{2 \mathrm{L}_{2}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right)^{(\mathrm{n}+1) / \mathrm{n}} \]

\[ \times\left[\left(\frac{\mathrm{R}_{2}\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{2}\right)}{2 \mathrm{L}_{2}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right)^{2} \frac{\mathrm{n}}{3 \mathrm{n}+1}+2 \unicode[Times]{x03C4}_{\mathrm{y}}\left(\frac{\mathrm{R}_{2}\left(\text{P}_ \mathcal a-\mathrm{P}_{2}\right)}{2 \mathrm{L}_{2}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right) \frac{\mathrm{n}}{2 \mathrm{n}+1}+\unicode[Times]{x03C4}_{\mathrm{y}}^{2} \frac{\mathrm{n}}{\mathrm{n}+1}\right] \]

\[ +\unicode[Times]{x03C0} \text{R}_{3}{ }^{3} k^{-1 / \text{n}}\left(\frac{\text{R}_{3}\left(\text{P}_ \mathcal a-\text{P}_{3}\right)}{2 \text{L}_{3}}\right)^{-3}\left(\frac{\text{R}_{3}\left(\text{P}_ \mathcal a-\text{P}_{3}\right)}{2 \text{L}_{3}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right)^{ \text{(n+1)} / \text{n}} \]

\[ \times\left[\left(\frac{\mathrm{R}_{3}\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{3}\right)}{2 \mathrm{L}_{3}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right)^{2} \frac{\mathrm{n}}{3 \mathrm{n}+1}+2 \unicode[Times]{x03C4}_{\mathrm{y}}\left(\frac{\mathrm{R}_{3}\left(\mathrm{P}_ \mathcal a-\mathrm{P}_{3}\right)}{2 \mathrm{L}_{3}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right) \frac{\mathrm{n}}{2 \mathrm{n}+1}+\unicode[Times]{x03C4}_{\mathrm{y}}^{2} \frac{\mathrm{n}}{\mathrm{n}+1}\right] \tag{3.15} \]

For Power Law flows, the relation corresponding to Equation 3.15, namely Equation 3.2, was solved by a “half-step” method to produce a numerical solution for the junction pressure P_a. Then, Equation 3.1 for general Q = {nπR³/(3n + 1)} [R(P_i – P_o)/(2kL)]^1/n was written for branch “1” in the form “Q₁ = {nR₁³/(3n + 1)} [R₁ (P_i – P_a)/(2kL₁)]^1/n ” and analytically solved to give the inlet pressure in Equation 3.6.

For Herschel-Bulkley fluids, the basic expression for flow rate is not so simple. Equation 3.15 is still solved as before in the Power Law case for the junction pressure P_a, but an analytical inlet pressure solution such as that in Equation 3.6 is not possible because of pressure nonlinearities. However, a numerical solution can be obtained if we write Equation 3.14 for branch “1,” that is,

\[ \frac{\text{Q}_{1}}{\unicode[Times]{x03C0} \text{R}^{3}}=k^{-1 / \text{n}}\left(\frac{\text{R}_{1}\left(\text{P}_{i}-\text{P}_ \mathcal a\right)}{2 \text{L}_{1}}\right)^{-3}\left(\frac{\text{R}_{1}\left(\text{P}_{i}-\text{P}_ \mathcal a\right)}{2 \text{L}_{1}}-\unicode[Times]{x03C4}_{y}\right)^{ \text{(n+1)} / \text{n}} \]

\[ \times\left[\left(\frac{\mathrm{R}_{1}\left(\mathrm{P}_{\mathrm{i}}-\mathrm{P}_ \mathcal a\right)}{2 \mathrm{L}_{1}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right)^{2} \frac{\mathrm{n}}{3 \mathrm{n}+1}+2 \unicode[Times]{x03C4}_{\mathrm{y}}\left(\frac{\mathrm{R}_{1}\left(\mathrm{P}_{\mathrm{i}}-\mathrm{P}_ \mathcal a\right)}{2 \mathrm{L}_{1}}-\unicode[Times]{x03C4}_{\mathrm{y}}\right) \frac{\mathrm{n}}{2 \mathrm{n}+1}+\unicode[Times]{x03C4}_{\mathrm{y}}^{2} \frac{\mathrm{n}}{\mathrm{n}+1}\right] \tag{3.16} \]

and solve for P_i also using a half-step approach. In summary, Q₁ is prescribed, while Q₂ and Q₃ are given by the first-two and second-two lines, respectively, on the right side of Equation 3.15. As before, the shear stresses are

\[ \unicode[Times]{x03C4}_1=\left(\text{P}_{i}-\text{P}_ \mathcal a\right) \frac{\text{R}_{1}}{2 \text{L}_{1}} \tag{3.17} \]

\[ \unicode[Times]{x03C4}_2=\left(\text{P}_ \mathcal a-\text{P}_{o, 2}\right) \frac{\text{R}_{2}}{2 \text{L}_{2}} \tag{3.18} \]

\[ \unicode[Times]{x03C4}_3=\left(\text{P}_ \mathcal a-\text{P}_{o, 3}\right) \frac{\text{R}_{3}}{2 \text{L}_{3}} \tag{3.19} \]

These follow from the global force balance “πR²ΔP = $\unicode[Times]{x03C4}$_w(2πRL),” noting that $\unicode[Times]{x03C4}$_y yield stress effects have already been considered in the calculations for pressure in Equations 3.15 and 3.16.

As before, we proceed with validations and re-run Examples 1-5 from our Power Law results, first for the zero yield stress case. Because the Power Law and Herschel-Bulkley equations and algorithms are completely different and solved by independently written programs, agreement with prior results, and with earlier Newtonian results, guarantees correctness. BIFURC-6 (not shown) provides the full suite of solutions for P_a, P_i, volume flow rate and viscous shear stress with yield effects (this solver also applies to Bingham Plastics for all k and $\unicode[Times]{x03C4}$_y values if we set n = 1).

3.2.2 BIFURC-6 Runs Assuming τ_y = 0 psi (Power Law limit)

C:\>bifurc-6

Example 1-5? Enter: 1

Solving for junction pressure PA ...

   N           PA                  QTEST           QCCPS           RATIO

369   0.200370E+01   0.9905E+02   0.1000E+03   0.9904621

370   0.200371E+01   0.9931E+02   0.1000E+03   0.9931461

371   0.200372E+01   0.9958E+02   0.1000E+03   0.9958305

372   0.200373E+01   0.9985E+02   0.1000E+03   0.9985147

Solving for Q1 inlet pressure P1 ...

   N           PINLET                  QTEST           QCCPS           RATIO

1477   0.201112E+01   0.9906E+02   0.1000E+03   0.9905897

1478   0.201113E+01   0.9913E+02   0.1000E+03   0.9912611

1479   0.201113E+01   0.9919E+02   0.1000E+03   0.9919319

.

.

1480   0.201114E+01   0.9926E+02   0.1000E+03   0.9926032

1490   0.201119E+01   0.9993E+02   0.1000E+03   0.9993138

1491   0.201119E+01   0.1000E+03   0.1000E+03   0.9999846

Additional calculated quantities ...

P1 = 2.0111935 (PSI)

Q1 = 100.00, Q2 = 50.06, Q3 = 50.06 (CC/S)

TAU1 = 0.00018659, TAU2 = 0.00009325, TAU3 = 0.00009325(PSI)

Stop - Program terminated.

The above results agree to the fourth decimal place with those in Figure 3 obtained using the exact Newtonian model and with those in Section 3.1 using the Power Law model. The results of Examples 2–5 here also agree with those in Section 3.1 to about four decimal places. This validates the “zero yield stress” limit of our solutions.

Example 1-5? Enter: 2

Solving for junction pressure PA ...

   N           PA                  QTEST           QCCPS           RATIO

.

1148   0.201150E+01   0.9966E+02   0.1000E+03   0.9966335

1149   0.201151E+01   0.9979E+02   0.1000E+03   0.9978742

1150   0.201152E+01   0.9991E+02   0.1000E+03   0.9991152

Solving for Q1 inlet pressure P1 ...

   N           PINLET                  QTEST           QCCPS           RATIO

3714   0.203011E+01   0.9903E+02   0.1000E+03   0.9902756

3715   0.203012E+01   0.9907E+02   0.1000E+03   0.9906565

.

.

3736   0.203022E+01   0.9987E+02   0.1000E+03   0.9986679

3737   0.203023E+01   0.9990E+02   0.1000E+03   0.9990498

3738   0.203023E+01   0.9994E+02   0.1000E+03   0.9994319

3739   0.203024E+01   0.9998E+02   0.1000E+03   0.9998142

Additional calculated quantities ...

P1 = 2.0302393 (PSI)

Q1 = 100.00, Q2 = 50.02, Q3 = 50.02 (CC/S)

TAU1 = 0.00046797, TAU2 = 0.00028802, TAU3 = 0.00028802(PSI)

Stop - Program terminated.

Example 1-5? Enter: 3

Solving for junction pressure PA ...

   N           PA                  QTEST           QCCPS           RATIO

.

12080   0.201152E+01   0.9997E+02   0.1000E+03   0.9997062

12081   0.201152E+01   0.9998E+02   0.1000E+03   0.9998245

12082   0.201152E+01   0.9999E+02   0.1000E+03   0.9999426

Solving for Q1 inlet pressure P1 ...

   N           PINLET                  QTEST           QCCPS           RATIO

.

39252   0.203024E+01   0.9999E+02   0.1000E+03   0.9998869

39253   0.203024E+01   0.9999E+02   0.1000E+03   0.9999232

39254   0.203024E+01   0.1000E+03   0.1000E+03   0.9999595

39255   0.203024E+01   0.1000E+03   0.1000E+03   0.9999962

Additional calculated quantities ...

P1 = 2.0302417 (PSI)

Q1 = 100.00, Q2 = 50.00, Q3 = 50.00 (CC/S)

TAU1 = 0.00046797, TAU2 = 0.00028807, TAU3 = 0.00028807(PSI)

Stop - Program terminated.

Example 1-5? Enter: 4

Solving for junction pressure PA ...

   N           PA                  QTEST           QCCPS           RATIO

.

19350   0.219376E+01   0.9991E+02   0.1000E+03   0.9991276

19360   0.219386E+01   0.9999E+02   0.1000E+03   0.9998652

19361   0.219387E+01   0.9999E+02   0.1000E+03   0.9999392

Solving for Q1 inlet pressure P1 ...

   N           PINLET                  QTEST           QCCPS           RATIO

.

62900 0.250880E+01   0.9998E+02   0.1000E+03   0.9998283

62905 0.250883E+01   0.9999E+02   0.1000E+03   0.9999418

62906 0.250883E+01   0.1000E+03   0.1000E+03   0.9999645

62907 0.250884E+01   0.1000E+03   0.1000E+03   0.9999871

Additional calculated quantities ...

P1 = 2.5088384 (PSI)

Q1 = 100.00, Q2 = 50.00, Q3 = 50.00 (CC/S)

TAU1 = 0.00787401, TAU2 = 0.00484695, TAU3 = 0.00484695(PSI)

Stop - Program terminated.

Example 1-5? Enter: 5

Solving for junction pressure PA ...

   N           PA                  QTEST           QCCPS           RATIO

.

8781   0.208793E+01   0.9994E+02   0.1000E+03   0.9993613

8782   0.208794E+01   0.9996E+02   0.1000E+03   0.9995888

8783   0.208795E+01   0.9998E+02   0.1000E+03   0.9998168

Solving for Q1 inlet pressure P1 ...

   N           PINLET                  QTEST           QCCPS           RATIO

.

24842 0.221233E+01   0.9997E+02   0.1000E+03   0.9997455

24843 0.221233E+01   0.9998E+02   0.1000E+03   0.9998258

24844 0.221234E+01   0.9999E+02   0.1000E+03   0.9999064

24845 0.221234E+01   0.1000E+03   0.1000E+03   0.9999870

Additional calculated quantities ...

P1 = 2.2123463 (PSI)

Q1 = 100.00, Q2 = 50.00, Q3 = 50.00 (CC/S)

TAU1 = 0.00310980, TAU2 = 0.00219886, TAU3 = 0.00219886(PSI)

Stop - Program terminated.

3.2.3 BIFURC-6 Runs Assuming τ_y= 0.00001 psi

Having validated our zero yield stress results, we turn to simulations with non-zero yield stress. We will re-run Examples 1-5 with a yield of 0.00001 psi, which is consistent with certain medical diagnostic examples, for example, Lee, Xue, Nam, Lim and Shin [5]. This larger value is not too much more than the values observed clinically but was chosen for numerical study purposes.

C:\>bifurc-6-yield-stress

Example 1-5? Enter: 1

Solving for junction pressure PA ...

   N           PA                  QTEST           QCCPS           RATIO

.

383   0.200424E+01   0.9923E+02   0.1000E+03   0.9923480

384   0.200425E+01   0.9950E+02   0.1000E+03   0.9950320

385   0.200426E+01   0.9977E+02   0.1000E+03   0.9977157

Solving for Q1 inlet pressure P1 ...

   N           PINLET                  QTEST           QCCPS           RATIO

.

1514   0.201224E+01   0.9976E+02   0.1000E+03   0.9975594

1515   0.201224E+01   0.9982E+02   0.1000E+03   0.9982303

1516   0.201225E+01   0.9989E+02   0.1000E+03   0.9989014

1517   0.201225E+01   0.9996E+02   0.1000E+03   0.9995728

Additional calculated quantities ...

P1 = 2.0122540 (PSI)

Q1 = 100.00, Q2 = 50.02, Q3 = 50.02 (CC/S)

TAU1 = 0.00019984, TAU2 = 0.00010651, TAU3 = 0.00010651(PSI)

Stop - Program terminated.

Example 1-5? Enter: 2

Solving for junction pressure PA ...

   N           PA                  QTEST           QCCPS           RATIO

.

1160   0.201202E+01   0.9971E+02   0.1000E+03   0.9971333

1161   0.201203E+01   0.9984E+02   0.1000E+03   0.9983738

1162   0.201204E+01   0.9996E+02   0.1000E+03   0.9996150

Solving for Q1 inlet pressure P1 ...

   N           PINLET                  QTEST           QCCPS           RATIO

.

3759   0.203126E+01   0.9986E+02   0.1000E+03   0.9985799

3760   0.203126E+01   0.9990E+02   0.1000E+03   0.9989616

3761   0.203127E+01   0.9993E+02   0.1000E+03   0.9993436

3762   0.203127E+01   0.9997E+02   0.1000E+03   0.9997258

Additional calculated quantities ...

P1 = 2.0312748 (PSI)

Q1 = 100.00, Q2 = 50.04, Q3 = 50.04 (CC/S)

TAU1 = 0.00048085, TAU2 = 0.00030102, TAU3 = 0.00030102(PSI)

Stop - Program terminated.

Example 1-5? Enter: 3

Solving for junction pressure PA ...

   N           PA                  QTEST           QCCPS           RATIO

.

12202   0.201204E+01   0.9997E+02   0.1000E+03   0.9997334

12203   0.201204E+01   0.9999E+02   0.1000E+03   0.9998513

12204   0.201204E+01   0.1000E+03   0.1000E+03   0.9999695

Solving for Q1 inlet pressure P1 ...

   N           PINLET                  QTEST           QCCPS           RATIO

.

39496   0.203127E+01   0.9999E+02   0.1000E+03   0.9998897

39497   0.203127E+01   0.9999E+02   0.1000E+03   0.9999259

39498   0.203127E+01   0.1000E+03   0.1000E+03   0.9999623

39499   0.203127E+01   0.1000E+03   0.1000E+03   0.9999988

Additional calculated quantities ...

P1 = 2.0312746 (PSI)

Q1 = 100.00, Q2 = 50.00, Q3 = 50.00 (CC/S)

TAU1 = 0.00048088, TAU2 = 0.00030098, TAU3 = 0.00030098(PSI)

Stop - Program terminated.

It is interesting to compare Example 1 (where the yield stress is zero) with Example 1 here (where the yield is 0.00001 psi). In the former case, P_a and P_i are 2.00373 and 2.01119 psi, while here, they are 2.00426 and 2.01225 psi. It is more meaningful to subtract the background “2 psi,” to compare 0.00373 with 0.00426 and 0.01119 with 0.01225. The percent increases are 14% and 9% respectively, which are large considering that $\unicode[Times]{x03C4}$_y = 0.00001 psi.

Example 1-5? Enter: 4

Solving for junction pressure PA ...

   N           PA                  QTEST           QCCPS           RATIO

.

19371   0.219437E+01   0.9998E+02   0.1000E+03   0.9998180

19372   0.219438E+01   0.9999E+02   0.1000E+03   0.9998918

19373   0.219439E+01   0.1000E+03   0.1000E+03   0.9999657

Solving for Q1 inlet pressure P1 ...

   N           PINLET                  QTEST           QCCPS           RATIO

.

62927 0.250986E+01   0.9999E+02   0.1000E+03   0.9999125

62928 0.250986E+01   0.9999E+02   0.1000E+03   0.9999353

62929 0.250987E+01   0.1000E+03   0.1000E+03   0.9999578

62930 0.250987E+01   0.1000E+03   0.1000E+03   0.9999806

Additional calculated quantities ...

P1 = 2.5098739 (PSI)

Q1 = 100.00, Q2 = 50.00, Q3 = 50.00 (CC/S)

TAU1 = 0.00788689, TAU2 = 0.00485995, TAU3 = 0.00485995(PSI)

Stop - Program terminated.

Example 1-5? Enter: 5

Solving for junction pressure PA ...

   N           PA                  QTEST           QCCPS           RATIO

.

8791   0.208843E+01   0.9994E+02   0.1000E+03   0.9993680

8792   0.208844E+01   0.9996E+02   0.1000E+03   0.9995956

8793   0.208845E+01   0.9998E+02   0.1000E+03   0.9998235

Solving for Q1 inlet pressure P1 ...

   N           PINLET                  QTEST           QCCPS           RATIO

.

24862 0.221333E+01   0.9997E+02   0.1000E+03   0.9997498

24863 0.221333E+01   0.9998E+02   0.1000E+03   0.9998301

24864 0.221334E+01   0.9999E+02   0.1000E+03   0.9999106

24865 0.221334E+01   0.1000E+03   0.1000E+03   0.9999914

Additional calculated quantities ...

P1 = 2.2133467 (PSI)

Q1 = 100.00, Q2 = 50.00, Q3 = 50.00 (CC/S)

TAU1 = 0.00312231, TAU2 = 0.00221136, TAU3 = 0.00221136(PSI)

Stop - Program terminated.

In closing, we note that the number of evaluations “N” has exceeded 60,000 in some runs. In others not reported here, they approach 200,000. This large number arises because we chose extremely small increments to obtain high accuracy solutions. Nonetheless, all cases required less than one second on Windows i5 machines. In our work, it was important to evaluate stability in the half-step approach, plus robustness and the ability to capture four decimal place accuracy in pressures and viscous stresses.

3.3 Newtonian and Herschel-Bulkley Examples

Blood flows in the aorta and larger arteries satisfy Newtonian models for which the viscosity μ is constant. The same viscosity value applies regardless of flow rate, vessel radius, position within the artery, or even changes in cross-sectional shape. This quantity is easily measured and its role is clearly understood. For instance, the Hagen-Poiseuille formula Q =πR⁴ (P_i – P_o) /(8 μL) illustrates how doubling μ will reduce the flow rate Q by half or how doubling the pressure gradient will double Q. These benefits, however, are lost in non-Newtonian flows. Three examples are considered.

For Example 1 for Newtonian flow, we show in Table 1 how our exact solutions provide clinical insights useful in studying flow blockage. We assume μ = 10 cp, Q₁ = 100 cc/s, R₁ = 1.52 cm, L₁ = 40 cm, R₂ = R₃ = 0.89 cm, L₂ = L₃ = 20 cm, and P₂ = P₃ = 2 psi. We have tabulated P₁, P_a,$\unicode[Times]{x03C4}$₁ and $\unicode[Times]{x03C4}$ ₂ (=$\unicode[Times]{x03C4}$₃) for simple clogs with decreasing flow areas (found straightforwardly by reducing radii successively by 10%) to model worsening blockages (other parameters are also easily considered). The results in Runs 8-10 show how severe consequences are expected at the higher blockage levels. (Of course, irregular blockages would be more realistic than concentric ones, but our purpose is to highlight the severity associated with large obstacles.)

Table 1 Example 1 – Newtonian results for aorta and iliac arteries.

As blood vessel diameters decrease, non-Newtonian effects become important. The behavior of the “apparent viscosity η,” or the proportionality factor between viscous shear stress and shear rate in “$\unicode[Times]{x03C4}$=ηdγ/dt,” is less clear. Unlike its Newtonian counterpart, η depends on flow rate, pressure drop, conduit size and shape, and so on. This means that blood viscosity will vary with local radius, pressure gradient, flow rate, geometric deformations and position in the cross-section. Thus, the η obtained from a specific lab viscometer shear rate may not be directly relevant to a diagnosis, so that the use of apparent viscosity for comparative assessments can be ambiguous.

However, the use of “n and k,” and “n, k and $\unicode[Times]{x03C4}$_y” for Power Law and Herschel-Bulkley fluids is unambiguous, because these constants represent true properties of the fluid. Together with Q₁, P₂, P₃ and the input geometry, parameters that will vary from case to case, the quantities P₁, P_a, $\unicode[Times]{x03C4}$₁, $\unicode[Times]{x03C4}$₂ and $\unicode[Times]{x03C4}$₃ particular to each scenario can be derived. In clinical work, it is more meaningful to study properties with dynamical importance rather than ill-defined viscosities with limited application. Preferred properties include pressure levels that affect bursting, pressure gradients related to flow efficiency, wall shear stresses that influence long-term wear and tear, and possibly ratios like Q₂/Q₃ that may affect overall body functions.

In Example 2 for Power Law fluids, we focus less on viscosity and more on dynamical properties as they depend on n and k. We assume Q₁ = 100 cc/s, R₁ = 1.52 cm, L₁ = 40 cm, R₂ = R₃ = 0.89 cm, L₂ = L₃ = 20 cm, P₂ = P₃ = 2 psi and $\unicode[Times]{x03C4}$_y = 0 psi. Typical results, based on published n and k data, are presented in Figure 10 (in the pressure plots, the background 2 psi is removed to highlight only dynamic effects). The plotting index L = 1 to 5 corresponds to n = 0.6 to 1.0, while that for M = 1 to 4 corresponds to k = 0.0000007325 to 0.000002930 (note that L = 5, M = 1 corresponds to a reference 5 cp Newtonian fluid). Herschel-Bulkley solutions for $\unicode[Times]{x03C4}$_y = 0.000001 psi showed very slight pressure and stress increases, so the corresponding figures are not plotted.

Figure 10 Example 2 – Power Law results for P₁, P_a, $\unicode[Times]{x03C4}$₁ and $\unicode[Times]{x03C4}$₂ (= $\unicode[Times]{x03C4}$₃).

Non-Newtonian effects, including yield stresses, are more pronounced in narrower blood vessels. In Example 3 for Herschel-Bulkley fluids, instead of repeating Example 2 with a nonzero $\unicode[Times]{x03C4}$_y value, we consider a system with smaller radii, both with and without yield stress. Comparative results for $\unicode[Times]{x03C4}$_y = 0 and a large value of $\unicode[Times]{x03C4}$_y = 0.00001 psi given in Figure 11 show relative numbers. The calculated wall shear stresses in all cases are within the ranges cited in published clinical data.

Figure 11 Example 3 – Herschel-Bulkley yield stress effects (tabulated results are given because plotted differences are visually imperceptible).

Power Law model limitations are well known in the literature. Note that the physical dimensions of “k” depend on the assumption for “n” – thus, for graphs like those in Figure 10, k values must not be compared when n values differ. In our work, the plots suggested another deficiency in the Power Law model. The formula Q = {nπR³/(3n + 1)} [R(P_i – P_o)/(2kL)]^1/n, reducing to Hagen-Poiseuille’s law with k =μ when n = 1, is well accepted in Newtonian flow. But suppose that Q is fixed, as is the case in blood flow studies. As n → 0, use of L’Hospital’s Rule leads to a questionable pressure gradient with (P_i – P_o)/L = 2k/R that is unrealistically independent of Q₁. This can be validated by direct calculations. Thus, caution is suggested in rate-constrained problems with small values of n – fortunately, n will typically exceed 0.7 in blood flow analyses. However, the same may not be true of some engineering calculations.

4. Discussion and Conclusions

While the subject of this paper, relating to non-Newtonian flows in bifurcated networks, appears to be one of fundamental importance, a detailed search through the civil, mechanical and process engineering literature found few analytically rigorous studies that were not empirically motivated. Our methods are based on sound mathematical and physical concepts. Because flows with nonlinear rheologies are each associated with their own special complicated details, we importantly described all the steps used in algebraic derivations and numerical methods. These were validated using exact solutions and cross-checks to guarantee that accurate solutions were produced. The authors hope that the body of formulas and procedures developed here will find good practical application in both engineering and biological sciences.

Acknowledgments

Gratitude is extended to the Reviewers for suggesting useful changes and clarifications. Also, colleagues at the Beijing International School, Beijing, China, are acknowledged for their interest, enthusiasm and constructive criticism.

Author Contributions

Jamie Chin developed the Newtonian models for pipe bifurcations, as well as the numerical analysis behind Power Law flows. Wilson Chin extended the methods to Herschel-Bulkley fluids, noting that while the overall strategies behind the use of “single pipe” solutions (in problems with multiple branches and for different rheologies) are similar, the implementations are quite different. The overall extension to conduits with arbitrary geometric clogs, also developed by Wilson, is described in detail in [13].

Funding

This work was internally supported by Stratamagnetic Software, LLC.

Competing Interests

The authors affirm that no competing interests exist. The research was entirely supported by internal Stratamagnetic Software, LLC funding. Portions of this paper were adapted from the authors’ forthcoming book Biofluids Modeling – Methods, Perspectives and Solutions, to appear in 2022, from scientific book publisher John Wiley & Sons.