A novel, FFT-based one-dimensional blood flow solution method for arterial network.

In the present work, we propose an FFT-based method for solving blood flow equations in an arterial network with variable properties and geometrical changes. An essential advantage of this approach is in correctly accounting for the vessel skin friction through the use of Womersley solution. To incorporate nonlinear effects, a novel approximation method is proposed to enable calculation of nonlinear corrections. Unlike similar methods available in the literature, the set of algebraic equations required for every harmonic is constructed automatically. The result is a generalized, robust and fast method to accurately capture the increasing pulse wave velocity downstream as well as steepening of the pulse front. The proposed method is shown to be appropriate for incorporating correct convection and diffusion coefficients. We show that the proposed method is fast and accurate and it can be an effective tool for 1D modelling of blood flow in human arterial networks.

Introduction

Mathematical and numerical modelling of blood flow in a human arterial network allows researchers to understand various flow related phenomena and disorders. It can help us to understand the genesis and progression of various diseases in the cardiovascular system, and it also provides us with a platform for developing methods for detecting such diseases. The arterial haemodynamics is influenced essentially by pulse wave propagation phenomena. At present time, one of the most effective and fruitful ways to understand wave phenomena in an arterial network is through one-dimensional (1D) flow modelling. Recent works on 1D blood circulation modelling have demonstrated its accuracy in predicting various flow quantities (Mynard and Nithiarasu 2008; Swillens et al. 2008; Alastruey et al. 2011; Gamilov et al. 2014; Sazonov et al. 2017).

The core of the 1D modelling is based upon numerical solution of nonlinear partial differential equations (PDE) derived, for example, in Formaggia et al. (2001), Sherwin et al. (2003). There are many numerical schemes/time integration methods available, including Finite Difference (FD) (Olufsen et al. 2000; Smith et al. 2002; Reymond et al. 2009; Saito et al. 2011), Finite Element (FE) (Formaggia et al. 2003; Sherwin et al. 2003; Mynard and Nithiarasu 2008; Alastruey et al. 2011), Finite Volume (FV) (Cavallini et al. 2008; Delestre and Lagrée 2013) and Discontinuous Galerkin (DG) methods (Matthys et al. 2007; Marchandise et al. 2009; Alastruey et al. 2011). Some of them are compared in the benchmark paper by Boileau et al. (2015). A rather fast scheme is proposed in Carson and Van Loon (2017) which is based on an Enhanced Trapezoidal rule method (Kroon et al. 2012). Many of these works appear to use slightly different mathematical models and constitutive relations as shown later in the present work.

Majority of the existing models, however, assume the convection velocity to be average over the cross section of an artery. Although this approximation is generally accepted, such approximations appear to give incorrect results, even for a steady-state flow. This error is pronounced in highly time-dependent flows, such as the pulsating flow experienced by an arterial network. In addition, space–time methods do not account correctly for the vessel skin friction of an arbitrary pulsating wave in an artery. The viscous effects should be based on the use of the Womersley solution (Womersley 1955) that is too complicated for the 1D space–time approach to handle. Other unsolved issues of 1D space–time approach include difficulties in implementing multi-elements Windkessel and/or non-reflective outlet boundary conditions, incorporate a robust viscoelastic vessel wall model, curvature of arteries and mass loss in smaller arteries and vessel walls.

An attempt to account for the vessel skin friction more accurately is made in the work by Bessems et al. (2007) where the authors consider differential velocity values near the wall region (a viscous layer) and the central core. Thus, the authors have introduced an additional degree of freedom. This approach can improve the accuracy of the computed skin friction to a certain extend but is less accurate than an approach based on the use of the Womersley solution. Another approach to account accurately for the skin friction in a space–time scheme is proposed in Reymond et al. (2009). The authors utilize a waveform obtained at a previous hear-beat cycle, perform the fast Fourier transform (FFT) for every element, build the Womersley solution by computing the Bessel functions of complex argument for every harmonic component, calculate contribution of every harmonic to the skin friction and perform the inverse Fourier transform. For every subsequent cycle, the accuracy of computed skin friction approaches to that of the Womersley solution. This method requires additional computations for every element.

An alternative to space–time method is proposed in Flores et al. (2016), which is based on linearization of the 1D equations and expanding the solution using the Fourier series. In this approach, the problem of wave propagation in an arterial network is solved analytically in the frequency domain, separately for every harmonic component. The pressure and flow rate waveforms are then calculated at any point of the network by numerically computing the inverse Fourier transform to the analytical solution obtained. In this approach, the skin friction can be accurately incorporated via the Womersley solution. Also, the viscoelastic properties of the vessel wall can be included without handling the complexities of the governing equations (Alastruey et al. 2011). As indicated in Flores et al. (2016), this approach allows one to rapidly investigate the role of individual physical properties of a cardiovascular system subjected to a pulsatile waveform. It is important to note that the Fourier transform and calculations carried out in frequency domain are the natural ways for dealing with the wave phenomena. Such important concepts as phase and group velocities can be employed in these methods to explain the distinctive features of wave the pulse propagation and reflection. The frequency domain approach has been successfully used in Sazonov et al. (2017) for developing a non-invasive, aortic aneurysm detection method using a waveform analysis.

Although the method proposed in Flores et al. (2016) is an excellent progress in terms of speed, it has many restrictions. To obtain an analytical solution using this method, every vessel (or its segment between two junctions) is approximated by a cylindrical pipe of constant cross section. Thus, the solution obtained is not general for realistic tapering vessels, found in arterial networks. Another restriction is that the linear algebraic equations have been derived manually for solving the linear problem and thus the method is not easy to employ on an arbitrary arterial network. Finally, the nonlinear effects are not considered by Flores et al. (2016). Therefore, despite its advantages, it is not competitive against the existing 1D modelling methods that use computational algorithms, especially in terms of robustness and accuracy.

In the present work, we propose the generalization of the FFT method by developing an effective and accurate procedure for solving the equations for a tapering vessel or vessel with abnormalities such as aneurysm and stenosis. Also, we account for nonlinear terms in the 1D equations through a nonlinear correction to the linear solution. In this way, we can capture the effect of increasing pulse wave velocity (PWV) downstream as well as steepening of the pulse front. In addition, in the implementation of the proposed method, the system of algebraic equations for every harmonic component is built automatically for any arbitrary arterial network. Finally, the nonlinear equations we employ contain the correct convection and diffusion coefficients.

This paper is organized into the following sections. In Sect. 2, the governing 1D arterial network equations are presented and the nonlinear and viscous terms in these equations are analyzed. Here, different constitutive relations are examined and compared. In Sect. 3, the proposed perturbation method is described, which reduces the nonlinear equations to a set of linear partial differential equations (PDE) which, in turn, are reduced to the ordinary differentials equations (ODE) by the Fourier transform method. An effective method for integration of the derived ODEs for an arbitrary tapering vessel is presented in Sect. 4 and generalized to rapid geometrical variations in vessels in the same section. In Sect. 5, the boundary condition required at inlet, outlet and vessel junctions are described. Here also the method for obtaining a linearized solution for an arbitrary arterial network is explained. In Sect. 6, a method for computing the second-order nonlinear corrections is presented for a single tapering vessel, for the boundary conditions and for full solution in an arterial network. Section 7 compares the present numerical results to those obtained by established 1D numerical methods and experimental results. In the final section, we outline the advantages and future prospects of the proposed approach. Some auxiliary material is presented in Appendix, including generalization of the Womersley solution for flow in a flexible pipe (“Appendix 1”).

Governing equations for 1D arterial network

The following system of partial differential equations (PDEs) for an arterial network can be derived from the incompressible Navier–Stokes equations (Formaggia et al. 2001; Sherwin et al. 2003)Here p(x, t) is the pressure; q(x, t) is the flow rate; A(x, t) is the lumen cross-sectional area; subscripts t and x denote partial derivatives with respect to time t and axial coordinate x, respectively; and are the blood density and kinematic viscosity, respectively; and are dimensionless coefficients defined and discussed in Secs. 2.1 and 2.2. Various forms of the constitutive relations are considered in Sect. 2.3.

Equations (1)–(3) form a basis for the 1D blood flow modelling. Note that these equations do not account for losses in the vessel wall due to its viscoelastic properties and also for blood flow rate losses in small lateral arterial vessels. In the following subsections we make some remarks on nonlinear and viscous terms in the governing PDEs.

The parameter

The coefficient is a factor in the nonlinear convection term in (2) and is defined aswhere u is the flow velocity profile; the integral is taken over the lumen cross section S. The coefficient is called the Coriolis parameter in Formaggia et al. (2001) and is also known as Boussinesq coefficient (Simakov and Kholodov 2009). Its value is assumed to be unity in most of existing models (Mynard and Nithiarasu 2008). This is actually valid for a flow with a uniform velocity profile , which has an infinite velocity gradient at the wall and, hence, the infinite value of the wall shear stress (WSS). Note that as shown in Sherwin et al. (2003), only for , (1)–(2) can be rewritten in terms of the mean velocity aswhich are mainly used in the numerical modelling works (eg., Mynard and Nithiarasu 2008; Gamilov et al. 2014). A steady-state flow in a cylindrical blood vessel of constant radius , i.e., Poiseuille flow, has the parabolic velocity profile of . This profile gives as shown in Vignon and Taylor (2004). In a pulsating flow, the instant velocity profile varies with time and so should the parameter (see Formaggia et al. 2003). An effective value of 1.1 for the parameter is estimated in Smith et al. (2002) and it is used, for example, in Alastruey et al. (2012a).

Friction coefficient

The friction coefficient for a cylindrical vessel with lumen radius of a can be calculated via the equationHere, subscript r denotes the partial derivative with respect to r, which is the polar coordinate in the lumen cross section S (). A uniform velocity profile gives as .

For the Poiseuille flow , the parameter takes a value of 4, which is used in most of the works (e.g. Mynard and Nithiarasu 2008). However, a pulsating flow is characterized by the high near-wall velocity gradient much higher than in the steady-state flow. Therefore an effective value of the coefficient, higher than four is used as well. For example, a value of 11 for is used in Alastruey et al. (2012a, b). This value is calculated in Smith et al. (2002) by fitting experimental data presented in Hunter (1972).

Analysis based on the Womersley solution (Womersley 1955) shows that the near-wall velocity gradient varies in time and can be higher or lower than that in the steady-state flow. Examples of time-variation of the and coefficients are shown in Fig. 1. They are computed using the 3D modelling of the flow in a Carotid artery described in Sazonov et al. (2011). The nonlinearity coefficient varies around a value of 4/3 and drops down to the value close to 1.1 only in vicinity the waveform peak. The instantaneous value of the friction coefficient has a larger range as shown in Fig. 1 and it reaches negative values under flow reversal conditions.

Fig. 1

Variation of the nonlinearity coefficient (red) and the friction coefficient (divided by 4) (blue) over a cardiac cycle. Normalized flow rate waveform is shown in black. These results are computed using the 3D modelling of the flow in a Carotid artery described in Sazonov et al. (2011)

From these results we can conclude that assuming coefficients and to be constant is very approximate. Assumption of an coefficient for an uniform velocity profile and the coefficient for the Poiseuille velocity profile in majority of existing studies should be changed to represent a more accurate pulsatile nature of the flow. A more rigorous approach can be based on convolution with functions and in the corresponding terms of the equations. These functions can be determined from the Womersley solution. However, this would complicate tremendously the space–time numerical schemes employed in the 1D modelling approach. Due to the reasons stated above, the nonlinear and viscous terms used in (2) by many of the space–time approaches are not very accurate. However, proposed method can easily incorporate different forms of convection and diffusion coefficients.

The constitutive relation

Now we consider the first term in (1). Applying the chain rule we can rewrite (1) in the formFrom here one can see that unless A(p) is a linear function, the first term is nonlinear with respect to p. Derivative can be referred to as the cross-sectional vessel compliance. The following dependence p(A) is proposed in Formaggia et al. (2001) and Sherwin et al. (2003)Here is the external/equilibrium pressure; is the lumen cross-sectional area at equilibrium state Sherwin et al. (2003); is a coefficient proportional to the vessel wall stiffness with h being the wall thickness and being the plate/shell analogue of the Young’s modulus and is the Poisson ratio. For an incompressible material and therefore . In Eq. (9), the parameter is assumed to be independent of A, which is valid if the vessel wall stiffness is assumed to be independent of strain.

Resolving relation (9) with respect to A we see that function A(p) is obviously nonlinear, i.e.,Some other dependencies of p(A) are considered below, and they are also nonlinear. Nevertheless, all of them can be expanded around the equilibrium state: :where prime denotes the derivative with respect to p and subscript 0 means that the value is calculated at the equilibrium state. In order to further analyse pressure–area relationship, lets introduce the pulse wave velocity (PWV) of a small perturbation to the equilibrium state aswhere is the lumen radius and subscript 0 denotes that the value is taken at the equilibrium state. Now, the cross-sectional vessel compliance can be represented using the following expansionwhere the ratiois a dimensionless form of pressure perturbation to the equilibrium state andare dimensionless nonlinear coefficients. They are calculated below for several types of constitutive law described in the literature.

Constant wall stiffness model Differentiating function (10) and using expansion (13) we obtainOlufsen’s model Model proposed by Olufsen (1999) (used in, eg, in Vignon and Taylor 2004) can be written in the form:This model givesi.e., threefold increase in nonlinearity coefficient compared to the constant stiffness model and additional nonzero coefficients .

Power law Both the above models are particular cases for the power law model proposed in Mynard et al. (2010), i.e.This relationship gives the constant stiffness model (10) for and Olufsen’s model (16) for . For the power law model,In Mynard and Smolich (2015) the following formula for the b parameter is proposed:where is the vessel collapse pressure. For example, for (indicated in Mynard and Smolich 2015), the reference pressure , and for a typical wave velocity in aorta, , we obtain , , , etc. Note that if b parameter obeys Eq. (20) or a similar law, then the nonlinearity parameters are different in different parts of the arterial system.

In the power law relationship, if then and the vessel wall becomes stiffer when the vessel expands during the pulse propagation. This dependence is indicated by works on arterial wall elastic properties (Holzapfel et al. 2000; Ogden and Saccomandi 2015). If then and the A(p) dependance becomes linear up to the third-order terms with respect to . There are other similar models in which the p(A) dependence is described in terms of elementary functions.

Armentano et al’s model One of the earliest models is proposed in Armentano et al. (1995). Here, the lumen diameter, , where and are constants. In terms of and the pressure-area relationship can be written as:which givesAs observed, the nonlinearity coefficients in the above equation are not constant. They can be very large for small arteries with high wave velocity.

Kholodov’s model An exponential/log dependence is proposed in Kholodov (2001) and employed in Gamilov et al. (2014), Simakov and Kholodov (2009), Vassilevski et al. (2015). For the dependence is exponential and is described bywhich givesHere, is negative and are nonzero.

It appears that not all studies use identical pressure-area relationships, and the nonlinear parameter varies from positive to negative values depending on the model used. These variations between studies indicate that the accuracy of the nonlinear part of the dependence p(A) is not well established and needs further investigations. Computational results of different studies can agree with each other only if nonlinear effects are small. A small nonlinear term in the A(p) dependence can be accounted through a correction to the linear problem. Thus, instead of the detailed study on dependance A(p) or p(A), one should focus only on nonlinear coefficients , especially , in large arteries, where pressure and area variations during the heart cycle are relatively large.

Note that integrating expansion (13), we obtain the useful relationshipThis indicates that the pressure and area amplitude parameters and are equivalent in the linear approximation.

Matching conditions at vessel junctions

At vessel junctions we enforce the law of mass conservation, i.e., the flow rate at the outlet of the parent vessel segment, , is equal to the sum of the inlet flow rates for the daughter vessel segments, . The mass conservation at vessel junctions may be written aswhere is the set of daughter vessel segments connected to parent segment s. Equation (25) is linear.

The continuity of the momentum flux requires that the total pressure should be continuous at vessel junctions (Sherwin et al. 2003), i.e.This is the matching condition used by majority of 1D arterial network models (Sherwin et al. 2003; Matthys et al. 2007; Mynard and Nithiarasu 2008; Marchandise et al. 2009; Alastruey et al. 2009) and it is nonlinear. Thus, the vessel junction matching condition is the third source of nonlinearity of the problem in addition to the convection term and pressure-area relationship. It is also important to note that some authors apply the continuity of the static pressure p, rather than the total pressure as the matching condition (Olufsen 1999; Reymond et al. 2009). A correction to the dynamic part of the pressure in the matching condition which depends on the angle between the parent and daughter vessels is discussed in Formaggia et al. (2003). The possibility to account for an additional pressure loss at junction is considered in Mynard and Valen-Sendstad (2015). This additional term, responsible for the pressure loss in the matching condition, is nonlinear with respect to velocity.

In summary, we can conclude that there are significant differences between studies in dealing with nonlinearity and often the nonlinearity is not satisfactorily represented. These variations are acceptable as long as the effect of nonlinearity is small. Since majority of the existing computational models assume limited nonlinearity in the arteries, we believe that a linearized perturbation method with nonlinear correction term is suitable for solving blood flow equations in a human arterial network.

Reduction of governing nonlinear PDEs to linear ODEs

Perturbation method

Multiplying (2) by and using the conservation mass equation (8) and also swapping the equations we can rewrite them in the formwhere is the dynamic viscosity. Substituting , into these equation and performing some manipulations we obtainNow retaining the linear and linearized terms on the left-hand side and substituting in virtu of (12), we, respectively, haveWe can rewrite Eqs. (27)–(28) in the matrix form:where is the linear operator, is the state vector, i.e., is a vector of nonlinear terms:If the nonlinear terms are small, Eq. (29) can be solved using subsequent approximation (perturbation) method , where satisfies the equationHere are the nonlinear correction terms:To derive (33) we have employed the following approximationsEquation (32) is linear for every m. It is homogeneous for . The right-hand side for depends on calculated at previous iterations: :

The Fourier transform method

The real heartbeat waveform is quasi-periodic, but the numerical simulations assume periodicity for convenience. Often simulations are repeated over several cardiac cycles before the results are used to avoid the effect of any initial conditions (Mynard and Nithiarasu 2008). Periodicity may typically be expressed asIn this case the Fourier series expansion is an appropriate method for solving the partial differential equations (PDE), i.e.,where i is the imaginary unity,is the frequency of nth harmonic component and and , respectively, represent the nth harmonic component (complex amplitudes of a harmonic wave at current x). They can be calculated using the integral over the cardiac period T, which is the direct Fourier transform of a periodic function, i.e.,Applying the Fourier transform (35) to Eq. (32), we obtain a system of ordinary differential equations (ODE):where is the nth harmonic component of the state vector . The linear operator contains only full derivatives, i.e.,The right-hand side of Eq. (36) becomesNow, instead of a set of PDEs, we have a set of independent linear ODEs for , where N is the number of harmonic components required. We do not need to solve these equations for negative n as a real results in , where denotes a complex conjugate.

Integration of the ODEs

Since we are working with the basic state values of A, q and c here, we omit subscript 0 in this section and in Sect. 6 as it can be confused with the zeroth harmonic component. In these sections we follow and and and regard them as independent of the wave amplitude. In the subsections below we consider the linear approximation, i.e. calculation of and the nonlinear corrections are considered in Sect. 6.

Zeroth harmonic component

In arteries of circular cross section, the flow is Poiseuille in nature and . Therefore, ODE (36) for takes the formwhich has the solutionwhere and are, respectively, the flow rate and the inlet pressure () averaged over a cardiac cycle, i.e.,

Nonzeroth harmonic components

In this section we omit subscript n in variables , and to make the equations easy to follow. For the linear approximation, we have the following ODEs:whereis a viscous factor which accounts for wave decays due to viscous losses, , and . The influence of viscosity on wave propagation and evaluation of is considered in Sect. 4.3. Explicit formulas for are derived in “Appendix 1” for a flexible cylindrical pipe.

Vessel of constant cross section and wall stiffness

The ODEs (40)–(41) do not admit an analytical solution in a closed form for a generic tapering vessel. Therefore, we first consider the simplest case of uniform vessel with constant parameters along the vessel axis.

Solution as a sum of travelling waves If A, c and are constant along the vessel, then we can write a general analytical solution in the linear approximation using the forward, , and backward, , harmonic travelling waves as particular solutions, i.e.,Here is the wavenumber; and are pressure amplitudes of, respectively, forward and backward harmonic waves; parameter is the characteristic vessel compliance that is an inverse value to the characteristic vessel impedance (Westerhof et al. 2005; Hametner et al. 2013)i.e. the pressure to flow rate ratio in the forward propagating wave along a uniform pipe. This parameter is frequency dependent in a pipe with viscous losses. This dependence is accounted in the factor . In solution form (43), the pressure and flow rate fields are fully determined by two numbers: and .

Effect of blood viscosity on wave propagation In the presence of viscosity, the wavenumber is complex . The real part of affects the wave speed (slowing it down) and increases the total phase incursion in the segment of length L. Its imaginary part, , indicates the wave decay during propagation. The wave decay per unit of length equals to .

Dependence of on the lumen radius a and the wave frequency is shown in Fig. 2(left) for calculated via Eq. (130) (see “Appendix 1”). One can see that its contribution to the wavenumber is essential for narrow vessels and at low frequencies. The wavenumber varies essentially along a tapering vessel at low frequencies. Therefore contribution of viscous losses looks to be more crucial for lower frequencies.

Fig. 2

Dependence of the imaginary part of the parameter (left) and the total phase incursion kL (right) on lumen radius a and frequency f (indicated above every curve in Hz). Here is taken

Due to viscous losses, the wave amplitude decays by a factor after propagation through the vessel. Therefore, since the total phase incursion is small for low frequencies, the relatively large imaginary part of the wavenumber does not cause the essential decay of the wave in a vessel. This is indicated in Fig. 2(right) where plots of are depicted for . One can see that the decay at the total length L is higher at a higher frequency.

Transmission matrix An alternative form of a general solution can be used if the vessel inlet conditions, P(0), Q(0), are given asHere is the transmission matrix, which is equal to the identity matrix at :In this approach, the pressure and flow rate fields are determined by two numbers P(0) and Q(0). Note that in Flores et al. (2016), two numbers: P(0) and P(L) are selected to determine the total field in the vessel. Below we show that the choice of P(0) and Q(0) has some advantages for numerical implementation.

Tapering vessel

ODEs (40) and (41) do not admit a closed from analytical solution for arbitrary functions A(x) and c(x). After analysing different possible approaches, we have selected the following two approaches.

Numerical integration of initial value problems ODEs (40) and (41) can be integrated numerically. To compute entries of transmission matrix , we should integrate ODEs (40) and (41) twice with the initial conditions, i.e.,The solution to the first initial value problem gives matrix elements and and the solution to the second initial value problem gives matrix elements and . This method is accurate and effective for the lowest harmonic components but looses its speed and accuracy for high-frequency harmonic components as the solution becomes fast oscillating.

Piecewise conic approximation Eliminating Q from ODEs (40) and (41), we obtain the second-order ODE with respect to P. After obtaining a solution, we can calculate Q asNote that wave velocity c varies along the vessel approximately proportional to , i.e. slower than the lumen radius. The variable varies fast only in a narrow and strongly tapering vessel. If we approximate c, and by their averaged values, we haveNow ODE (48) admits an analytical solution in a number of cases. We consider here the case of a conic vessel segment: , then we can write the general solution in the form of forward and backward waves in which the pressure amplitude grows toward the narrower edge of the vessel, i.e.,whereThe transmission matrix for this vessel iswhere , , .

A tapering vessel can be approximated by a sequence of truncated conic elements with values and k, approximated by their average value over every element. Let be edge points of the cones. Here, is the number of conic elements.

To accelerate the computation, the inlet and outlet values of every ith element can be used to calculate the corresponding mean values:where subscript i indicates that the value is calculated at point . Using (52) we can write the transmission matrix of every element asHere,where is the lumen radius at . The transmission matrix of total segment can be calculated through the product of transmission matrices for every element as:Accuracy of the piecewise conic approximation To study the accuracy of the piecewise conic approximation, we compute the reflection and transmission coefficients of a tapering vessel, located between two cylindrical pipes with matched areas and wave velocities. The Right Carotid segment has been selected from the arterial network described in Mynard and Nithiarasu (2008) having simultaneously the largest length and highest lumen radius gradient: the length is , the inlet and outlet diameters, respectively, are and . The vessel is approximated by a cone. The wave velocity dependence on lumen radius is taken to be approximated by Eq. (102) proposed in Blanco et al. (2015).

Three basic element lengths are selected and for the conic approximation. These lengths are adjusted to the vessel length by the rule: , , where means the closest integer equal or greater than x. Therefore the element lengths are adjusted to values and .

The reflection and transmission coefficients computed in the conic approximation are compared with the same coefficients , computed by the numerical integration with the accuracy of . Results of comparison are presented in Fig. 3. Observe that the accuracy is very good, especially for the first five harmonic components containing more than 90% of the pulse wave energy.

Fig. 3

Accuracy of the reflection coefficient (left) and transmission coefficient versus frequency f for different element lengths h (indicated above every curve)

Boundary conditions

Consider a network containing vessels and vessel segments with one inlet segment and a number of terminating segments. In every segment the axial coordinate x is referenced from its inlet. The flow parameters and in the sth segment can be calculated by the methods described above if the inlet flow parameters and are known.

Introduce the notations and , where . Thus are inlet amplitudes for every segment of the network. To compute the flow in the entire network, we have first to calculate unknowns , for every harmonic component. Hence, we have to derive equations with respect to these unknowns from the inlet, junction and outlet boundary conditions.

In Flores et al. (2016) the following unknowns are taken: . This is suitable if the vessels are approximated by uniform pipes without nonlinear correction. In stead of dealing with and , we can actively use the transmission matrices to calculate solutions for tapering vessels. Also, this approach simplifies computation of nonlinear corrections.

Inlet boundary condition Let the periodic pressure waveform is set at the inlet of the network . We regard as a full pressure taken from direct measurements (i.e. sum of forward and the backward waves). Then the boundary conditions iswhere are calculated via the Fourier transform (35).

If is a forward propagating waveform (as it used in Mynard and Nithiarasu (2008) as an inlet boundary condition), then we can rewrite solution (50) in the formwhere and are inlet parameters of the first element in the first segment. Now eliminating the amplitude of the backward wave , we obtain an equation with respect to unknowns and . Substituting and , we havewhere is the characteristic impedance of very first element of the first segment.

Analogously, if the periodic flow rate waveform is given at the inlet and it is the full waveform then the inlet boundary condition iswhere are calculated via the Fourier transform (35).

If is the only forward propagating waveform then for the nonzeroth harmonics, we can derive in the similar mannerMatching conditions at junctions Consider the sth segment of length connected to daughter segments at its outlet. Remind that and are inlet flow parameters for the sth segment for the nth harmonic component. Now, the outlet flow parameters can be written in the linear approximation through the transmission matrix for the sth segment asHereinafter without argument means . For the zeroth harmonic relations (61) are different (see solution (39))but they can be written in the form of (61) if we define the transmission matrix for , i.e.,If we denote the set of daughter segments associated with a parent segment s as , the continuity of the total pressure and the mass conservation at junction can, respectively, be written asSubstituting relations (64) into conditions (61), we obtain linear homogeneous equations asMatching conditions at junctions with merging arteries In some junctions two parent arteries can merge to form a single daughter artery. For example, the left and right vertebral arteries merge to form the basilar artery which is a common daughter segment for the vertebral arteries. There are few similar backward bifurcations in the cerebral arterial system. Similar situation can occur in an arterial system with a bypass. In a general case, for the node with parent segments and daughter segments, we have the following set of equationsThe number of equations for the junction always equals the nodal index of the junction: . In the particular case of a junction with occurring among the cerebral arteries gives us three equationsThese equations are valid for zero harmonic too if the transmission matrix is taken in the form of (63).

Boundary condition at the outlets of terminating segments In general, the Windkessel model is employed at the outlet of a terminating vessel (Alastruey et al. 2012a; Boileau et al. 2015; Carson and Van Loon 2017), which is analogous to an electric circuit containing resistors and capacitors. In the present approach proposed, for every harmonic component with frequency , the impedance load for a terminating segment s can be defined by single complex number rather than by an ODE as in the traditional numerical schemes, i.e.,For example, in the standard three-element Windkessel model, containing two resistors and and a capacitor C, the impedance equals .

The non-reflecting boundary condition can be set in a natural way for any nonzeroth harmonic component. In the proposed approach, the outlet impedance should be equal to the ratio of the pressure and flow rate in the harmonic forward propagating wave, i.e.,If the terminating segment is cylindrical, then the outlet impedance simply should be equal to the characteristic impedance for that segment: . If the terminating segment is approximated by a set of cones (see Sect. 4.4), then as it is seen from (50), the P / Q ratio in the forward propagating wave is . Therefore the non-reflecting outlet impedance should beHere superscript ‘out’ denotes the value at the outlet of terminating segment s. Expressing and through the transmission matrix we obtain the linear equation for and asValue in condition (71) has sense of characteristic impedance of a conic vessel with constant wave velocity. The load with such impedance will absorb all forward propagating waves.

Zeroth harmonic component The inlet/junction/outlet conditions should be considered separately for the zeroth harmonic component. The component for the inlet segment is set by Eq. (56). At junctions, we have the matching conditions (65). As for the outlet conditions at terminating vessel segments, the outlet resistance can be set at all such segments as:where is the pressure at which flow to the microcirculation ceases (Alastruey et al. 2012b).

Alternative way is to apply Murray’s law (Murray 1926; Sherman 1981), splitting the flow rates at a junction proportional to a power of lumen radius, i.e.,where . Some studies show that a value of is optimal for laminar flow, for turbulent flow and is a good choice for arterial blood flow (see Olufsen 1999 and references there in). If additional information is known about splitting the flows at some junctions for particular networks, it also can be utilized instead of (74) at those junctions.

Solution to the linearized problem Conditions (58) [or (60)], along with conditions (65) and (72) form a closed systems of linear algebraic equations with respect to and for all . Equations (60), (65) and (73) [or (74)] form a closed system for .

At this point, we have closed systems of linear equations with respect to and for and for all frequencies accounted: . Recollect that for a negative n we have the same solutions but complex conjugate. Then for monitored site of a selected vessel segment s, the waveform can be computed through the transmission matrix and the inverse Fourier transform, (34),The Fast Fourier Transform (FFT) method is very useful in this case. Examples of the solutions are presented in Sect. 7.5

Nonlinear corrections

Nonlinear correction to solutions Let inlet parameters and for the sth vessel segment are found for all n (subscript s is omitted in this section). The next stage is to substitute them into relation (33) in order to calculate the right-hand sides for equations at the second iteration. We should mention here that the most expensive part is to compute the parameter in relation (33).

Also observe that in Fig. 1 the variation of is relatively small around the value 4/3 in contrast to the larger variation of the parameter. Also, note that the accuracy in the computation of the correction to the main approximation can be lower compared to the main approximation. Therefore from practical viewpoint, the easiest way is to approximate by a constant value, say, 4/3. Accurate computation of the parameter and study of its influence on the waveform will be addressed in the future works.

Approximating by a constant we can rewrite relation (33), after omitting subscript 0 (indicating the basic state) and performing the differentiation of the product, asEstimations show that the third addend in the first row accounting for combined correction of nonlinearity and viscosity is relatively small and can be omitted. Utilizing properties of the Fourier transform and also Eqs. (40) and (41) we can calculate the derivatives asAfter computation of for the segment, we have to apply the Fourier transform (38) to relation (76) and integrate the inhomogeneous ODEs to givewith zeroth initial conditions at the segment inlet: and . The zeroth initial condition here needs explanation. Solution to linear inhomogeneous ODEs (80) and (81) with general initial conditions can be represented as a sum of solution to homogeneous ODEs (40) and (41) with general initial conditions (which gives solution ) and the solution to inhomogeneous ODEs with zeroth initial condition (this addend just gives solution ).

ODEs (80) and (81) can be integrated numerically, for example, using Runge–Kutta methods, but we propose here a less accurate but faster method of integration based on piecewise conic approximation. Recall again that a lower accuracy is acceptable for correction to the main approximation. Solution to the inhomogeneous ODEs with zeroth initial conditions can be represented through the quadratures asSubscript n is omitted here for brevity. To calculate entries of the matrix, we use the conic approximation described in Sect. 4.4. To calculate the quadratures, we employ the Trapezoidal rule using cones’ edge points as points of discretization. Now, the flow parameters at the outlet of the jth element are given aswhereHere , where is given by formula (54). In the case of , solutions to ODEs (80) and (81) are given asAgain the quadratures can be computed by the Trapezoidal rule asNonlinear corrections to the junction matching conditions With the account for nonlinear corrections, the outlet flow parameters take the following formwhere and satisfy inhomogeneous ODEs (80)–(81) with zero initial conditions. Now, we can invoke the continuity of total pressure at a junction. Equation (85) can now be rewritten aswhereare nonlinear terms accounted only in the second iteration after linear values, and , have been computed. Values are computed through the inverse and direct Fourier transforms:where . Here and are, respectively, the inlet and outlet area of the sth segment, and . Now, substituting (87), we obtain linear inhomogeneous equations asAnalogous equation can be written for a more complicated arterial junction in which more than one parent artery are connected.

Solving the nonlinear problem The second iteration of the network problem equations is different from the first iteration as inhomogeneous equations (90) instead of (65) are employed in the second iteration. If Murray’s law, (74), is used, then the corrected outlet flow rate also should be accounted. The total system of algebraic equations with respect to and is still linear, but most of the equations are inhomogeneous. After solving all systems of equations, the corrected waveforms can be computed in all necessary sites of a flow network.

Comparison with other numerical 1D schemes for an arterial network

Effect of accurate viscous term

First, we consider the effect of the more accurate account of viscous decay of the propagating pulse. For this purpose we compute propagation of a single Gaussian-shaped pulse along a uniform pipe having length of , lumen radius of and , considered first in Alastruey et al. (2012a) and then used for comparison of various numerical 1D schemes in Boileau et al. (2015). The flow rate at the inlet is set as /s. Non-reflection boundary condition is set at the outlet. The standard blood properties are used, i.e., and /s. The velocity profile at inlet is taken as  Alastruey et al. (2012a), Boileau et al. (2015) with that corresponds to .

The velocity reaches at a peak value of 1 cm/s which is much smaller than the wave velocity. Therefore nonlinear effects are small and we can apply the linearized approach described in Sects. 4 and 5. In our approach, we have to deal with periodic pulses. Thus, we select a period that is significantly greater than the time of pulse propagation along the pipe, . The results of computation with the frequency independent are plotted using black lines in Fig. 4 (left). One can compare these results with that of Alastruey et al. (2012a), Boileau et al. (2015) and confirm that pulse peak is identical to the published results.

Fig. 4

Left: propagation of a Gaussian-shape pulse along a 10 m pipe of the 2 cm diameter with . Black: ; Red: Womersley’s profile-based decay. Dashed line indicates the theoretical decay (93). Right: propagation of a realistic waveform along the same pipe at for (black solid), for (black dashed), and for Womersley’s profile-based decay (red)

Using the approach developed here, this decay can be evaluated analytically. The wave number that accounts for viscosity isFor the case under consideration that is small in comparison with frequency of all harmonic components except the zeroth one. Taking into account the fact that is small and expanding givesi.e. all harmonic components with frequency much higher than are decaying at the same rate asThus, if such harmonic components mainly constitute the pulse, then the pulse will propagate, preserving its shape and decaying exponentially with the coefficient . This coefficient is equal to for the problem considered here. The narrow pulse considered here is rich in high-frequency harmonics for which decay given by (93) is accurate. The analytical exponential decay of the pulse peak is plotted by a dashed line in Fig. 4.

When correctly accounted for viscous term, i.e. when the parameter is calculated accurately based on Womersley’s solution (see Flores et al. 2016 and Sect. A), change in shape and amplitude of the propagating pulse as shown in Fig. 4 (left) is observed. One can see here that the pulse is decaying noticeably faster, its shape is modified with formation of pulse tail, and the pulse peak propagates slightly slower than in the case of constant . Such big difference in pulse decay is obtained as we have considered a very narrow pulse instead of a realistic waveform, observed in arteries. This narrow pulse is rich in high-frequency harmonics with value greater than 11. In a realistic waveform, first five harmonic components contain 95% of the pulse energy, whereas in the narrow Gaussian pulse used in Alastruey et al. (2012a), Boileau et al. (2015), the number of harmonic components containing 95% of the pulse energy is 32. For a realistic waveform computed with the constant gamma and based on Womersley profiles give approximately the same decay of the pulse peak as seen in Fig. 4 (right). Nevertheless there is some difference in the pulse shape and propagation speed between the constant approximation and Womersley’s profile-based treatment of the viscous term. A large difference is observed in the narrow negative peak in Fig. 4 (right) as it is formed by higher frequencies.

For the sake of comparison, a waveform propagation computed with the (the value often taken in 1D arterial network modelling) is plotted using the dashed line in Fig. 4 (right). One can see that underestimates the pulse peak decay in the pipe with . This indicates that in a large artery like aorta is a reasonably effective viscous parameter, averaged over accurate values of the main harmonics constituting a typical waveform. At the same time, in a narrower vessel, the effective constant should be smaller than 11 as the Womersley’s profile-based values for the main harmonics are smaller (see Fig. 13).

Fig. 13

Dependence of real part (red), imaginary part (blue) and absolute value (black) of the parameter on the Womersley number . Dashed lines indicate an approximation (131)

Nonlinear effects in a uniform pipe

Comparison with the analytical solution First, we estimate the accuracy of the numerical integration (84) against the analytical solution for a uniform pipe and inviscid flow. Equation (32) may be written in the following form:Here the pressure and flow rate are expressed through dimensionless variables , . Equation (94) admits an analytical solution, for example, for the following linear wavewhere is an dimensionless amplitude factor: ; parameter determines the ratio of the constant and oscillating parts. Solutions to PDEs (94) satisfying boundary condition are:whereThe dimensionless nonlinear correction parameters along the 10-m pipe problem (see previous subsection) for the linear waves described by Eq. (95) are presented in Fig. 5. One can see that exact and numerically obtained curves are indistinguishable. Analysis shows that relative error of the numerical computation is 0.06% for an element length of .

Fig. 5

Dimensionless nonlinear corrections for the pressure (left) and the flowrate (right) for the waveform, (95), with . Solid lines: the exact solution, (96) and (97); dashed darker lines: results of the numerical calculations, (84) with

Solution (96)–(97) helps us to compare contribution of different nonlinear parameters in a waveform. The growing part of the solution contains the following factorThe convection term, , mainly contributes to and for the Poiseuille profile and the term in Eq. (27) is unity in relation (98). The value of in the constant stiffness constitutive law is 0.5. Thus, contribution of nonlinearity, related to the constitutive law, is the small in this case. Nevertheless the constitutive law (20) can give large negative value of . In such situation, contribution of the vessel wall to the nonlinear corrections can be much higher and it can additionally increase the growth of nonlinear corrections.

Accuracy for a tapering vessel

Here, we consider the simplest network which comprises a 1-m-long pipe with uniform cross section, connected to a tapering segment and the outlet of the tapering part is connected to another 1-m-long pipe. Such long inlet pipe is chosen to have forward and reflected pulses clearly separated in time (see Fig. 6). The dimensions of the tapering segment correspond to the Right Carotid segment used in the arterial network described in Mynard and Nithiarasu (2008) and used in Sect. 4.4. Its inlet and outlet diameters, respectively, are 6.75 mm and 3.5 mm, and its length is 9.4 cm. The shape of the tapering segment here is approximated using a truncated cone. The diameters of the first and third pipe are equal, respectively, to the inlet and outlet diameters of the tapering segment. The dependance of PWV on lumen diameter is taken from Blanco et al. (2015). The non-reflecting boundary condition (69)–(70) is set at the outlet of the network. A flow rate with a bell-like shape is set at the inlet of the network, i.e.,where the period, , and the pulse duration, .

Fig. 6

Pressure (left) and flow rate (right) waveforms computed at the inlet, 1, midpoint, 2, outlet, 3, of the first pipe and outlet of the conic segment, 4. The pressure waveform in cite 4 not presented as it is almost identical to that of at site 3. Coloured curves are computed using the method described in Carson and Van Loon (2017), dashed black curves are computed using the proposed method

The pressure and flow rate waveforms are computed at four sites: 1—the inlet, 2—midpoint and 3—outlet of the first pipe and 4—the outlet of the tapering segment. The waveforms computed by the numerical scheme described in Carson and Van Loon (2017) are shown in Fig. 6 using coloured lines and waveforms computed by the proposed method are shown using dashed black lines. The inlet flow amplitude is taken sufficiently small with /s to ensure that the process is linear and is employed in both computations. One can see that agreement between proposed method and numerical computations is excellent and this proves the accuracy of the proposed method, at least for the linear process.

We now discuss how the waveforms depend on the inlet peak flow amplitude, . A part of the waveforms in site 3 (inlet of the tapering segment), in the vicinity of the main peak are shown in Fig. 7. They are shown in a normalized form with and , which are proportional to the amplitude parameter, . Therefore, dividing by we obtain normalized waveforms having approximately the same amplitude distorted only by the nonlinear factors. One can see that the pressure peaks are shifted forward with the increase in amplitudes. The waveforms computed by the proposed method behaves similarly to that of the numerical computation, but their shift is smaller when the amplitude approaches /s.

Fig. 7

Part of the normalized pressure (left) and flow rate (right) waveforms computed at site 3 (inlet of the tapering segment) in the vicinity of the main peak for different peak flow rate values (indicated in the legend). The solid coloured lines are computed using the method in Carson and Van Loon (2017), and the dashed lines are computed using the proposed method

In Fig. 8 one can see the evaluation of nonlinear corrections and with the growth of the inlet flow amplitude. As the second-order correction are proportional to the square of the amplitude, it is natural to plot them in the normalized form , . In the proposed method, if we account only the second-order corrections in equation (32), the normalized waveform is amplitude independent. They are indicated by the black dashed line in Fig. 8. To compute the nonlinear corrections using standard 1D methods, in particular, using the method in Carson and Van Loon (2017), we carry out the following procedure. First, we compute the waveform for an extremely low amplitude. We take /s, compute waveforms and and use it as a reference linear solution. Then, we calculate the waveforms p and q at a realistic amplitude and calculate the nonlinear correction by the formula,Here , where is the mean pressure. The results of the computation are plotted in Fig. 8 in the normalized form (solid lines). The distortion of this corrections with the growth of amplitude indicates influence of the higher-order corrections in Eq. (32).

Fig. 8

Normalized nonlinear corrections for the pressure (left) and flow rate (right) waveforms computed at site 3 for different flow rate peak values of the inlet flow. The colours are the same as in Fig. 7. Normalized waveforms and computed using the proposed method is independent of the amplitude and indicated by black dashed curve

Note that the nonlinear distortion of the pulse is high in this network due to a long () pipe connected to a tapering segment. The consideration based on Riemann equation, describing formation of a shock wave confirms this. The nonlinear distortion is accumulated during its propagation along a long uniform pipe as the higher harmonics generated by the nonlinear effect propagate with the same speed along a uniform pipe. Any change in geometry like tapering, branching, PWV variation, etc. will cause discrepancy between the phase velocity of the different harmonics that will prevent accumulation of nonlinear distortion. Therefore, in the real arterial system such distortion should be smaller.

Comparison with experimental results

In the work Sazonov et al. (2017), an in vitro setup is described for studying pulse reflection from synthetic aneurysms of different diameters. Details of the setup, parameters, inlet flow rate waveforms and measured pressure waveforms are presented in Sazonov et al. (2017). A comparison of the measured pressure waveforms and waveforms computed by the proposed method is depicted in Fig. 9 by black and red lines, respectively. The waveforms are compared for the site from the inlet of the main pipe, i.e. at the midpoint of the 50-cm-long pipe preceding the segment with an aneurysm. The comparison is carried out for all aneurysm diameters, , considered in Sazonov et al. (2017), i.e., 24 mm, 34 mm, 44 mm and 50 mm.

Fig. 9

Comparison of measured (black) and computed (red) pressure waveforms. The experimental setup is described in Sazonov et al. (2017)

All the pipes in the experimental setup, except one with the aneurysm have a constant diameter. Therefore, the calculation of the transmission matrix is straightforward using (45). The element size for nonlinear corrections is taken equal to 2 cm. For a 14 cm segment with an aneurysm, the element size is taken equal to 0.5 cm, both for the linear problem and for nonlinear corrections. The viscosity factor is taken frequency dependent and is computed by the approximation, (130). The nonlinearity coefficients are: , . In the experiments described in Sazonov et al. (2017), single pulses are used because of the strong reflection from the outlet of the pipe network. As in the proposed method we have to deal with periodic signals, the period is taken large (4 s), and non-reflecting boundary condition is set at the outlet of the pipe network. The small discrepancy between computed and measured waveforms in Fig. 9 is caused by the uncertainty of mechanical parameters of some elements of the setup and uncertainty of the inlet boundary conditions.

Application of the method to arterial networks

In most arterial network models, the exact variation of the blood vessel parameters is not specified as subject-specific information is very difficult to obtain (see for eg. Mynard and Nithiarasu 2008). In many cases, we have to interpolate the vessel parameter variations along the length of the vessel. For example, the lumen area A and the wave speed c [or the parameter in relation (9)] require interpolation along the length. Often A(x) is linearly interpolated along the vessel length (Mynard and Nithiarasu 2008) and the same approach is followed for parameter if its inlet and outlet values are given. In some works, the authors approximate only A(x) and apply one of the known approximations for c(a) (recall is the lumen radius). Such approximations are mainly employed in 1D modelling, which can be written in the following forms:where the lumen radius a should be in centimetres, and the obtained PWV c is in m/s. Approximation (101) is proposed in Olufsen (1999) and its parameters for systemic arteries are listed in Mynard and Smolich (2015). Approximation (102) is proposed in Blanco et al. (2015) and approximation (103) is proposed in Reymond et al. (2009). Here, 20a equals to the lumen diameter D expressed in mm. These approximations for c(a) are plotted in Fig. 10(left).

Fig. 10

Left: The vs lumen diameter D approximations. Blue: Olufsen (1999), Mynard et al. (2010), Mynard and Smolich (2015); Red: Reymond et al. (2009); Green: Blanco et al. (2015). Right: Different types of vessel shape approximation. Blue: A(x) is linear. Red: a(x) is linear. Green: A, a are exponential

Different types of lumen radius variation, adapted by different authors, are depicted in Fig. 10(right). They are: the area is linearly interpolated as in Mynard and Nithiarasu (2008) (blue), the radius is linearly interpolated as shown in red and exponential interpolation of area is shown in green. One can see that the linear interpolation of the area produces a convex and non-realistic shape. Thus, other types of interpolation of A(x) are better approximations. Nevertheless, for a typical artery, the difference between inlet and outlet radii is much smaller and discrepancies between the interpolations are hardly important.

For assessing the accuracy of the proposed method, we compare our results against the numerical modelling results of the human arterial network using a well established method, described in Mynard and Nithiarasu (2008). Here we use the arterial network described in Mynard and Smolich (2015), which contains 107 blood vessel segments as shown in Fig. 11. The dependence of the PWV on lumen radius is calculated using Eq. (102). The flow rate is imposed at the inlet of the first segment as in Blanco et al. (2015); Boileau et al. (2015) with the heartbeat period of . All the tapering segments are approximated by truncated cones. We have set constant values and as in the computational model Mynard and Nithiarasu (2008). The constitutive relation is taken in the form of relation (9) and thus the nonlinearity parameter . The three-element Windkessel (lumped) model is employed as terminal boundary conditions.

Fig. 11

Arterial network model described in Mynard and Smolich (2015) and used here for simulation: main arteries (a), cerebral arteries (b), coronary arteries (c)

The comparison of results between the present and numerical computations is shown in Fig. 12. Observe that the proposed method gives results very close to that computed by the method used in Mynard and Nithiarasu (2008). There are some small discrepancies near the peak and in the decaying part of the pulse wave, which can be attributed to the computational errors. The numerical model of Mynard and Nithiarasu (2008) needs about 300 s to compute the first 3 heartbeat cycles. The proposed method needs only 6 s to compute the correct waveforms, including the nonlinear corrections. The model proposed in Carson and Van Loon (2017) needs 50 s to compute waveforms in the same networks. Overall, the proposed method is fast and as accurate as the numerical models mentioned.

Fig. 12

Pressure wave form (left) and flow rate waveform (right). Computed waveforms at the beginning of aortic arc (green), at the beginning of abdominal aorta (blue) and at the midpoint of the right carotid (red). Waveforms computed by the model in Mynard and Nithiarasu (2008) are shown using solid lines and those computed by the proposed method are shown by dashed lines

Discussion and conclusions

Many of the currently used numerical models for arterial flow do not account correctly for the viscous friction, especially for a pulsating flow. It is obvious that several different values of viscosity coefficient are employed by currently used models. The nonlinear convection term used in the mathematical models in the past is not rigorous for a pulsating flow. In addition, a variety of constitutive laws are employed for vessel walls, some of which resulted in nonlinear terms in the governing equations. All these variations in mathematical model suggest that there is no consensus among researchers on how to deal with nonlinear effects when solving for flow in arterial networks. Since the nonlinearity is often small, all the treatments used in the literature appear to work without too much disagreement in certain flow regimes (low amplitude waves). However, we have proposed a generalization of the nonlinear terms through corrections that may be applicable for all flow regimes.

A novel method, based on the perturbation technique and fast Fourier transform (FFT) in solving the 1D blood flow equations, has been proposed in the present work. The proposed method makes the FFT competitive against traditional space–time numerical schemes in terms of both robustness and speed. In contrast to the FFT approach described in Flores et al. (2016), the proposed method can be applied to an arbitrary arterial network, containing tapering vessels, vessels with stenosis and aneurysms, and to a high amplitude waveform in which the nonlinear effects are relevant. As demonstrated by the results, the proposed method is faster than competing methods and it is accurate. It accounts for viscous effects more accurately than any existing space–time methods and more importantly the viscous coefficient, , is automatically calculated for different flows and physical conditions. The proposed method simplifies boundary conditions required at the terminal vessels. Thus, we believe that this method can be an alternative and potentially a more effective tool for 1D modelling of blood flow in arterial networks.

Although the proposed method is a substantial improvement to the existing methods, it requires further development in the following area. Further attention is required to deal with viscoelastic effects, blood mass loss due to smaller branches, porous nature of arteries and application to clinical environment.