## Abstract

Dynamic analysis of large-scale flow networks is made difficult by the large system of differential-algebraic equations resulting from its modelling. To simplify analysis, the mathematical model must be sufficiently reduced in complexity. For self-similar tree networks, this reduction can be made using the network’s structure in way that can allow simple, analytical solutions. For very large, but finite, networks, analytical solutions are more difficult to obtain. In the infinite limit, however, analysis is sometimes greatly simplified. It is shown that approximating large finite networks as infinite not only simplifies the analysis, but also provides an excellent approximate solution.

## 1. Introduction

Potential-driven transport networks, such as mammalian circulatory and bronchial systems, heating, ventilating and air-conditioning systems, city plumbing networks, and power distribution grids, frequently appear in many branches of science and engineering. While the nature of these systems varies, the modelling of each network produces a system of differential-algebraic equations (DAEs) governing the flow. For small-scale systems with only a few branches and nodes, there are several techniques to solve the DAE system. But large-scale systems, or systems with a large number of nodes and branches, result in large systems of DAEs that can be expensive to solve.

There are many techniques that can be used for steady-state analysis of small-scale transport networks. The Hardy Cross method (Cross 1936) of analysing fluid piping networks and Kirchoff’s circuit laws for analysing electrical circuits are two such techniques. But dynamic analysis of even a small network is much more difficult and requires the solution of the DAE system (Mayes *et al.* in press). For large-scale transport networks even steady-state analysis can become complicated owing to the extremely large number of equations in the resulting mathematical model. In the same manner, as the scale and complexity of a transport network increase, time-dependent solutions can be increasingly difficult to obtain and the problem can become effectively intractable because of its size.

DAE systems are often classified by their index and while different definitions exist, all are a measure of the difficulties one can expect in calculating a solution. For large, high-index DAE systems, such as those that arise from transport networks, there are methods to reduce the index (Hedengren & Edgar 2005). At best the result of such methods is a large system of ordinary differential equations (ODEs)Q2 that must still be reduced or simplified. Oftentimes this reduction is in the form of assumptions or approximations, such as assuming one-dimensional flow or using linear approximations when modelling. Other methods of reduction work to reduce the size of the mathematical model. For example, a complex circuit can be reduced to a much simpler macromodel of the original circuit using Thévenin’s and Norton’s equivalent circuits. This reduces the number of equations in a mathematical model and can greatly reduce the computational expense of obtaining time-dependent solutions.

Similarly, some mathematical models can be reduced by taking advantage of physical structure or self-similarity in the original network to create a simpler network model. Nakagawa & Sorimachi (1992) do this in the case of an infinite resistor–capacitor circuit with a physically repeated pattern to create a simpler, finite macromodel circuit, allowing the dynamic solution of a previously intractable mathematical problem. Similar approaches have been used with other infinite electrical networks, including networks with differing forms of similarity such as ladders, grids and rings (Zemanian 1991; Srinivasan 1992; Mavromatis 1995; Thompson 1997; Kelly & McGough 2009). In each of these cases, models are reduced by taking advantage of physical structure and the infinite nature of the system.

All of these methods of reduction and simplification seek to reduce complexity in such a way that solutions can be obtained. In most cases, a system is reduced by approximating it as something smaller in scale or simpler in nature. But reduction of a network can also be accomplished by approximating it as a network even larger in scale if that approximation simplifies it in a way that allows solutions to be obtained more easily. In the present work, large-scale systems are reduced by approximating them by their infinite extensions which paradoxically allows for simpler solutions. Potential-driven flows in bifurcating tree networks are considered as examples to show how the system model can be reduced from a very large or infinite set of DAEs to a single ODE.

## 2. Trees

Bifurcating tree networks, also referred to as perfect binary trees and shown in figure 1, regularly occur in both natural and artificial engineering systems. Examples include biological transport networks, river basin drainages, viscoelasticity models and microchannel electronic cooling systems (Tarboton *et al.* 1988; Heymans & Bauwens 1994; Cross 1997; Cieplak *et al.* 1998; Chen & Cheng 2002; Pence 2002; Alharbi *et al.* 2003; Dokoumetzidis & Macheras 2003; Masters 2004; Gabrys *et al.* 2006; Wang *et al.* 2006; Escher *et al.* 2009). Its appearance in so many natural and man-made systems is no coincidence; Bejan (Dan & Bejan 1998; Bejan 2000) has shown the bifurcating tree network geometry to be the optimized result of point-to-volume flow systems.

This network begins with a bifurcation into two branches, which make up the first generation. Each branch in the first generation then splits into two new branches in the following generation. This repeats itself for each successive generation. A network of *N* generations then results in 2^{N+1}−2 total branches and 2^{N}−2 branching points. Each branch can be represented uniquely with a pair of integers (*i*, *j*), where 1≤*i*≤*N* is the generation number the branch is located in, and 1≤*j*≤*M*_{i}=2^{i} is the branch number within that generation.

The system studied is assumed to have three basic properties: flow (of fluid, heat, energy, etc.) occurs through the network; the transfer is driven by a potential difference across the system; and the transferred quantity is conserved at each of the branching points. With these assumptions, a transfer equation can be written for every branch of the form
2.1where *q*_{i,j}(*t*) is the transfer through branch (*i*,*j*), Δ*p*_{i,j}(*t*) is the potential difference across the branch (*i*,*j*) and *L*_{i,j} is the operator relating the two. The operators *L*_{i,j} could each be of any form, the only restriction made is that they be linear and, if necessary, restricted by appropriate initial conditions such that a unique inverse exists. For the operator *L* and the equation *Lx*=*y*, the solution *x* is given by
2.2where *L*^{−1} is the inverse of *L*, *α*_{n} is a constant and *ϕ*_{n} belongs to the null space of the operator *L*, i.e. *Lϕ*_{n}=0. For many engineering systems, either the kernel of the operator is trivial or that initial conditions can be chosen such that . In either case, the rightmost term in equation (2.2) can be disregarded. With these conditions, the inverse operator *L*^{−1} is unique.

In addition, a conservation equation can be written for each branching point within the network of the form
2.3where *P* is the number of branches leaving each branching point and *w*_{k} is a weighting factor.

For any bifurcating tree network this results in a set of 3×2^{N}−4 equations. In its most general form, this DAE system can be written as
2.4where *q* and Δ*p* are the vectors of flow rate and potential difference, respectively. Because some of the equations in this coupled system are algebraic, it must be noted that the matrix *A* is guaranteed not to be full rank and is thus not invertible. This is, by definition, a differential-algebraic system. Additionally, the matrix *E*, while relatively sparse, is not banded. As the size of the network grows, the resulting system of DAEs increases exponentially in size and complexity. For large-scale systems, with *N* very large, reduction then becomes not just useful, but necessary to obtain solutions. The goal is to approximate this large or infinite set of DAEs by a single ODE
2.5where *q*(*t*) is the total transfer through or across the network, Δ*p*(*t*) is the potential difference across the network and *L*_{N} is the approximate operator relating the potential difference and induced transfer. This will provide a reduced mathematical model to allow for dynamic solutions of even very large systems.

### (a) Reduction in the time domain

In an effort to reduce the mathematical model of a transport network, the goal is then to find the operator *L*_{N} in equation (2.5) which is equivalent to the complete tree composed of the operators *L*_{i,j} that describe the flow through each branch. The system of equations can, however, be simplified by noting the self-similar structure present in both the network and the resulting equations, and using that observation to reduce the size of the equation set and eliminate unnecessary unknowns. An equation of the form shown in equation (2.1) exists for every branch within the network. These equations can be added along any possible path from network inlet to outlet to eliminate the intermediate potential differences. For example, consider the 2-generation network shown in figure 2. There are 2^{N} (with *N*=2) possible paths through this network. The six transfer equations for this system (one for each branch) are
2.6a
2.6b
2.6c
2.6d
2.6e
2.6fwhere *q*_{i,j}, Δ*p*_{i,j}, and *p*_{i,j} are functions of time as before. Additionally, assuming unit weights the conservation equations for the two nodes are
2.7aand
2.7bFinally, the total flow, *q*, through the simplified network is given by
2.8and *L*_{2} is the operator describing the behaviour of the simplified 2-generation tree in
2.9

By combining the transfer equations along the four unique paths from inlet to outlet, the interior potentials, *p*_{1,1} and *p*_{1,2}, are eliminated to yield four new equations
2.10a
2.10b
2.10c
2.10dEquations (2.10*a*–*d*) can then be manipulated to solve for the *q*_{2,j}’s by moving the leftmost term to the right-hand side and then operating on both sides by , which is the inverse operator that when applied to *L*_{2,j} produces the identity operator, i.e. . The *q*_{2,j}’s are then
2.11a
2.11b
2.11c
2.11dUsing equations (2.7*a*,*b*), these equations combine to give
2.12aand
2.12bThe process is then repeated to solve for and eliminate the *q*_{1,j}’s by first applying the inverse operators and to equations (2.12*a*,*b*), respectively, to give
2.13aand
2.13bCollectingQ1 the *q*_{1,j}’s on the left and operating with and , respectively, results in
2.14aand
2.14bThis can then be substituted into equation (2.8), resulting in the single equation
2.15Rewriting equation (2.15) in the form of equation (2.9), the system operator for a bifurcating network with *N*=2 generations can be given as
2.16

This same process can be repeated for a tree network with *N* generations. The 2^{N+1}−2 branch equations can be summed along all of the 2^{N} possible unique paths through the network to eliminate all of the intermediate potentials. The conservation equations can then be used to eliminate the unknown flow rates within the network, beginning with the *N*th generation and progressively working back to the first, just as was done before. For a large *N*, this process results in the reduction of a very large set of equations to a single equation in terms of the overall flow rate and potential difference. For an *N*-generational potential driven tree system, *L*_{N} can be written as
2.17

### (b) Reduction in the frequency domain

The usefulness of equation (2.17) depends on how easily *L*_{N} can be evaluated. For algebraic operators, computation of the system operator is simple as the inverse operators are merely the algebraic inverse of the operator. But in the case of integral or differential operators, evaluation of equation (2.17) can be complicated even when it is assumed that all the inverse operators are unique. As many dynamic analyses of engineering systems are concerned with differential operations, this severely limits the utility of equation (2.17). However, while evaluation of *L*_{N} in the case of differential or integral operators may be difficult in the time domain, it can be calculated relatively easily in the frequency domain by use of the Laplace transform.

Consider the same two-generation system as shown before in figure 2. If it is assumed that each of the operators *L*_{i,j} are differential, then the set of DAEs (2.6*a*–*f*)–(2.9) results. Taking the Laplace transform of this set of DAEs and following the same procedure as before, the system operator can be written as
2.18where the notation ℒ_{i,j}=ℒ[*L*_{i,j}] denotes the Laplace transform of the operator *L*_{i,j}.

Furthermore, for differential or integral operators *L*_{i,j}, the Laplace transform of the operator is a simple polynomial in the frequency domain. The inverse operator then is simply the reciprocal of ℒ_{i,j}. Taking advantage of this, ℒ_{N} can be written in the form of a continued fraction,
2.19With ℒ_{N} expressed as a continued fraction a variety of techniques can be used to find either an exact or approximate convergent of the operator in the Laplace domain. Integrating back to the time domain then results in a simplified mathematical model consisting of a single ODE for the complete system.

For smaller systems, ℒ_{N} can be calculated for any ℒ_{i,j}’s. For large *N*, ℒ_{N} can either be calculated in the same way or be approximated as . For certain forms of similarity in the operators, this can result in very clean and simple expressions. Equations (2.18) and (2.19) in their current form are in the most general form possible for a bifurcating tree network where each of the operators *L*_{i,j} can be different and unrelated to one another. While all bifurcating tree networks have a structural self-similarity to them, there are additional forms of similarity, both within a generation and between successive generations of the network, which can allow additional simplification.

The first form of similarity considered is that within each generation. If the operators within the *i*th generation are identical then the system is said to be symmetric, otherwise it is asymmetric. The second form of similarity is between the operators of successive generations. If the operators present in the *i*th generation are dependent on either the generation number *i* or the operators in the preceding generation, then the system is said to be generation-dependent, otherwise it is generation-independent.

While the expressions seen in equations (2.18) and (2.19) are valid only for bifurcating tree networks, they can be generalized to *P*-furcating networks with
2.20and
2.21

Though the previous analysis is focused on trees, the concept is applicable to a much more general class of problems where modelling results in a very large set of DAEs. For an infinite tree, the problem becomes one of reducing and solving an infinite number of equations which could be much simpler than solving a large, but finite, equation set. The physical structure of the tree provides an easy method of spotting a pattern or similarity within the equation set which allows it to be reduced or simplified. Patterns or similarity modes can also be found in other large or infinite systems of equations that would allow reduction to a single equation regardless of whether the equation set is the result of a potential driven network, chemical reaction or mechanical system.

## 3. Examples and results

Any bifurcating tree network results in a system of DAEs. For extremely large networks, finding dynamic solutions of the resulting system of DAEs can often be computationally expensive or even impractical. But by taking advantage of structure to reduce the mathematical model, simple approximations and even exact analytical solutions can be found.

### (a) Piping networks

One such example is that of a piping network arranged in a bifurcating tree. Consider an *N*-generational bifurcating tree network of circular pipes transporting fluid owing to a potential difference across it. If the diameter and length of the pipe are taken to vary with generation number *i* according to *D**_{i,j}=*D**/*β*^{i−1} and ℓ*_{i,j}=ℓ*/*β*^{i−1}, where *D**=*D**_{1,j} and ℓ*=ℓ*_{1,j}, where * denotes a dimensional quantity, then appropriate values of *β* can be chosen to ensure that all physical quantities, such as length and volume, remain finite even for an infinite system (Franco *et al.* 2006).

Using a one-dimensional flow model for laminar pipe flow, the dimensional momentum equation for a branch is
3.1where *t** is time, *q** the flow rate, Δ*p** is the pressure drop across a pipe, *D** is the diameter, ℓ* the length of pipe, and *ν** and *ρ**are the kinematic viscosity and density of the fluid, respectively. Non-dimensionalizing by *q*_{i,j}=*q**_{i,j}/(*ν***D**^{2}), and *t*=32*ν**/*D**, and substituting for and , equation (3.1) becomes
3.2In the form of *L*_{i,j}*q*_{i,j}=Δ*p*_{i,j}, the generational-dependent operator can be written as *L*_{i,j}=*β*^{i−1}*d*/*dt*+*β*^{3(i−1)}, and in the Laplace domain, ℒ_{i,j}=*β*^{i−1}*s*+*β*^{3(i−1)}. Using equation (2.19) and solving for the convergent to the continued fraction in an *N*th-generation network,
3.3
3.4which is a simple polynomial in the frequency variable *s* that is dependent on *N*.

Transforming back from the Laplace domain an exact analytical solution for the behaviour of the entire network can be found,
3.5In this example, a large set of DAEs is reduced to a simple ODE, and an exact solution for a network of any size *N* can be found easily. But it is often the case that for very large *N*, finding a simple analytical expression such as equation (3.4) is not trivial or may not even be possible. Even if an analytical expression can be found, transforming back to the time domain may be challenging. While simple and exact solutions may not be possible, approximate solutions can be easily found. From equation (3.5), as *N* tends to infinity, solutions converge to the infinite network case. In this case, for 1<*β*<2^{1/3} (the same restrictions on *β* can be shown to be required for all network quantities to remain finite) the solution converges to
3.6

Figure 3 shows the step response of bifurcating piping networks for different values of *N* and the error between those responses and that of the infinite system. In this case, the error, *E*, is defined as
3.7This definition is chosen to best represent the difference between the dynamic responses of the finite and infinite solutions. It is obvious from figure 3 that as the size of the system increases, solutions converge to the infinite system response and that even relatively small systems can be well approximated by an infinite network. In this case, approximating the network as infinite provides a very good approximate solution even for networks of only *N*=8 or *N*=12 generations.

### (b) A fractance device

A typical RLCQ2 circuit is composed of capacitors, resistors and inductances which have impedances *Z*∝(*jω*)^{m}, where *m* is −1,0 and +1, respectively. A fractance is a device which does not fall into any of the three common cases, but instead has an impedance *Z*∝(*jω*)^{α}, where *α* is not an integer. Nakagawa & Sorimachi (1992) use an infinite tree of resistors and capacitors of the form of figure 1, with *L*_{i,j} replaced by resistors for *j*=odd and capacitors for *j*=even, to produce such a device. For a branch with a resistor, the voltage/current relationship is given by *Ri*(*t*)=*v*(*t*), where *R* is the resistance of the resistor, *i*(*t*) and *v*(*t*) are the current and voltage, respectively. The branches with capacitors are governed by
3.8where *C* is the capacitance of the capacitor and *i*(*t*), *v*(*t*) are again, the time-dependent current and voltage, respectively. In the notation of equation (2.1), *q*=*i*, Δ*p*=Δ*v*. Taking the operators to be *L*_{i,j=odd}=*R*(⋅) for the resistor branches and for the branches with capacitors, their frequency domain counterparts are given by ℒ_{i,j=odd}=*R* and ℒ_{i,j=even}=1/(*Cs*). Using these operators and finding the convergent of the infinite continued fraction (2.19) gives . The impedance as a function of *ω* is found by replacing the Laplace variable *s* with *jω*, giving *Z*=(*R*/*C*)^{1/2}(*jω*)^{−1/2}, which is indeed a fractance device. Converting back from the Laplace domain,
3.9which gives *v*(*t*) as exactly the 1/2th integral of *i*(*t*). Equation (3.9) can also be rewritten to show *i*(*t*) as the fractional order 1/2 derivative of the voltage.

Finding the infinite network solution in this case is simple. More challenging is to find solutions for large, but finite, networks. In this case, a simple and easily invertible expression for the convergent of the continued fraction in equation (2.19) is not available. For smaller networks () the convergent can be found and numerically inverted, although it is increasingly difficult to do so as the size of the network increases. For networks larger than *N*≈10, even the numerical Laplace inversion of equation (2.19) becomes prohibitively expensive. Figure 4 shows the voltage response to a step input in current for several small networks, as well as the step-response of the infinite case. In this case, the step-response of the infinite network is unbounded owing to its infinite resistance. As steady state is never reached, any error calculation would be somewhat arbitrary and as such is not included. But it can clearly be seen that as the size of the network increases, solutions approach that of the infinite network solution. For short times, small networks again provide an excellent approximation.

The fractance device problem becomes more difficult if the generational dependence of the pipe flow example is added. Assuming that the resistance and capacitance in generation *i* is given as *R*_{i}=*β*^{i}*R* and *C*_{i}=*C*/*β*^{i}, where *β*≤1/2 to ensure that the resistance of the infinite network remains finite, then in the frequency domain the operators are ℒ_{a}=*β*^{i}*R*, and ℒ_{b}=*β*^{i}/(*Cs*). For the infinite network, substituting into equation (2.19) gives a quadratic equation in ,
3.10

resulting in two possible solutions. Numerically inverting back to the time domain, both solutions can be found. In this case, the solution resulting from the negative root can be seen to be spurious as the solution does not behave as would be intuitively expected.

Figure 5 shows the step-response of the generationally dependent fractance device using the positive root of . As before, solving the infinite network case is simpler than solving for large finite networks. Assuming large networks to be infinite, however, again provides an acceptable approximate solution. Figure 5 shows the step-response of *N*=1, 2, 6, 8 and generation networks, and the error between the infinite and finite networks with increasing generation number *N*. The error is computed in the same manner as equation (3.7).

## 4. Conclusions

Transport networks are common in all branches of engineering, in all shapes and sizes, but dynamic analysis is often complicated by the complexity in the resulting mathematical model. Even for relatively small networks, dynamic solutions can be difficult to obtain owing to the size of the ensuing system of coupled DAEs. For very large networks, solutions can become effectively intractable making model simplification or reduction a necessity.

In the case of bifurcating tree networks, the self-similar structure of the system and its resulting mathematical model allows for massive reduction. It can be expressed in the frequency domain in terms of a continued fraction. And by finding the convergent of the fraction the large or even infinite system of DAEs from the original mathematical model can be reduced to a single ODE.

As was shown in the included examples, it may only take a relatively small number of generations before the response rate approaches that of the infinite network. In light of these examples, the question should then be asked why a very large system should be approximated as infinite when even approximating the system with as few as *N*=6 generations could yield an acceptable solution. This method—reduction of complexity by a reduction in size—is indeed the approach taken most often when modelling large systems. But it should be noted that the systems shown in this paper are exceptional cases chosen precisely because of the fact that they allow for simple solutions of smaller networks. Generally speaking, this is not the case. Consider that even an *N*=6 generation results in an order-2 DAE system of 188 coupled equations. For more complicated system operators, this is not trivial at all. And in these cases, reduction of complexity by approximating the system as infinite then has the advantage of producing both a more tractable and more accurate solution.

While the development in this paper focused on linear systems with a perfect binary tree structure, it should be noted that these requirements are not so restrictive as to make it inapplicable to real-world systems. A similar approach can be taken to analyse systems with other self-similar or nearly self-similar structures. While the requirement that the systems be linear is a strong one, it is still helpful for large scale nonlinear systems. In these cases, linearization is just another step in the reduction of a large-scale, nonlinear set of DAEs to a simpler mathematical model.

## Acknowledgements

J.M. acknowledges the support of the late Mr D. R. Dorini as well as the US Department of Energy and the Centre for Applied Mathematics at the University of Notre Dame.

- Received February 12, 2011.
- Accepted April 6, 2011.

- This journal is © 2011 The Royal Society