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), . 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) 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; 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 A = A(p) 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 as where 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 u ≡ū , 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 , only for = 1 , (1)-(2) can be rewritten in terms of the mean velocity ū = q∕A as which 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 a = √ A∕ , i.e., Poiseuille flow, has the parabolic velocity profile of u(r) = 2ū (1 − r 2 ∕a 2 ) . This profile gives = 4∕3 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 equation Here, subscript r denotes the partial derivative with respect to r, which is the polar coordinate in the lumen cross section S ( r ∈ [0, a] ). A uniform velocity profile gives = ∞ as u r (a) = ∞.
For the Poiseuille flow u r (a) = 4ū∕a , 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 u r (a) 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. 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 (t) and (t) 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 form From here one can see that unless A(p) is a linear function, the first term is nonlinear with respect to p. Derivative dA∕dp ≡ A � (p) can be referred to as the cross-sectional vessel compliance. The following dependence p(A) is proposed in Formaggia et al. (2001) and  Here p 0 is the external/equilibrium pressure; A 0 is the lumen cross-sectional area at equilibrium state (p, q) = (p 0 , 0)Sherwin et al. (2003) to the vessel wall stiffness E ′ h with h being the wall thickness and E � = E∕(1 − 2 ) being the plate/shell analogue of the Young's modulus and is the Poisson ratio. For an incompressible material = 0.5 and therefore E � = 4 3 E . In Eq. (9), the parameter is assumed to be independent of A, which is valid if the vessel wall stiffness hE ′ 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: p = p − p 0 : 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 c 0 as where a = √ A∕ 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 expansion where the ratio is a dimensionless form of pressure perturbation to the equilibrium state and are 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 obtain Olufsen's model Model proposed by Olufsen (1999) (used in, eg, in Vignon andTaylor 2004) can be written in the form: This model gives i.e., threefold increase in nonlinearity coefficient 1 compared to the constant stiffness model and additional nonzero coefficients 2 , 3 , ….
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)  where p collapse is the vessel collapse pressure. For example, for p collapse = −10 mmHg (indicated in Mynard and Smolich 2015), the reference pressure p 0 = 80 mmHg , and for a typical wave velocity in aorta, c 0 = 5 m/s , we obtain b = 4.33 , 1 = −1.17 , 2 = 1.94 , etc. Note that if b parameter obeys Eq. (20) or a similar law, then the nonlinearity parameters 1 , 2 , … are different in different parts of the arterial system.
In the power law relationship, if b > 1 then 1 < 1∕2 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 b = 2 then 1 = 0 and the A(p) dependance becomes linear up to the third-order terms with respect to p . 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, D = + ln p , where and are constants. In terms of p 0 , A 0 and c 0 the pressure-area relationship can be written as:

which gives
As 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), Kholodov (2009), Vassilevski et al. (2015). For A ≥ A 0 the dependence is exponential and is described by which gives Here, 1 is negative and 2 , 3 , … are nonzero.
It appears that not all studies use identical pressure-area relationships, and the nonlinear parameter 1 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 n , especially 1 , in large arteries, where pressure and area variations during the heart cycle are relatively large.
Note that integrating expansion (13), we obtain the useful relationship This indicates that the pressure and area amplitude parameters p = p∕ c 2 0 and Â = A∕A 0 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, q s , is equal to the sum of the inlet flow rates for the daughter vessel segments, q d . The mass conservation at vessel junctions may be written as where 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 p + 1 2ū 2 should be continuous at vessel junctions , i.e. This is the matching condition used by majority of 1D arterial network models 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.

Perturbation method
Multiplying (2) by ∕A 0 and using the conservation mass equation (8) and also swapping the equations we can rewrite them in the form where = is the dynamic viscosity. Substituting into these equation and performing some manipulations we obtain Now retaining the linear and linearized terms on the lefthand side and substituting A � 0 = A 0 ∕( c 2 0 ) in virtu of (12), we, respectively, have We can rewrite Eqs. (27)-(28) in the matrix form: where L is the linear operator, u is the state vector, i.e., n is a vector of nonlinear terms: If the nonlinear terms are small, Eq. (29) can be solved using subsequent approximation (perturbation) method Here n (m) are the nonlinear correction terms: To derive (33) we have employed the following approximations Equation (32) is linear for every m. It is homogeneous for m = 1 . The right-hand side for m > 1 depends on p (m � ) , q (m � ) calculated at previous iterations: m ′ < m:

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 as In 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 P n (x) and Q n (x) , 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 U (m) n = [P (m) n , Q (m) n ] T is the nth harmonic component of the state vector u (m) . The linear operator contains only full derivatives, i.e., The right-hand side of Eq. (36) becomes Now, instead of a set of PDEs, we have a set of independent linear ODEs for n = 0, 1, 2, … , N , where N is the number of harmonic components required. We do not need to solve these equations for negative n as a real u (m) n * , 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 c = c 0 and A = A 0 and q = q 0 and regard them as independent of the wave amplitude. In the subsections below we consider the linear approximation, i.e. calculation of P (1) n , Q (1) n 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 = 4 . Therefore, ODE (36) for n = 0 takes the form which has the solution where q and p(0) are, respectively, the flow rate and the inlet pressure ( x = 0 ) averaged over a cardiac cycle, i.e.,

Nonzeroth harmonic components
In this section we omit subscript n in variables P (1) n , Q (1) n and n to make the equations easy to follow. For the linear approximation, we have the following ODEs: where is a viscous factor which accounts for wave decays due to viscous losses, 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, Here k = ∕c is the wavenumber; C f and C b are pressure amplitudes of, respectively, forward and backward harmonic waves; parameter Ỹ = A∕( c ) 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: Effect of blood viscosity on wave propagation In the presence of viscosity, the wavenumber is complex k = ∕c . The real part of affects the wave speed (slowing it down) and increases the total phase incursion Re (kL) in the segment of length L. Its imaginary part, Im ≤ 0 , indicates the wave decay during propagation. The wave decay per unit of length equals to exp{−Imk}.
Dependence of −Im( ) on the lumen radius a and the wave frequency f = ∕2 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.
(43) Due to viscous losses, the wave amplitude decays by a factor exp{−ImkL} 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 −ImkL are depicted for L = 10 cm . 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 as Here T(x) is the transmission matrix, which is equal to the identity matrix I at x = 0: 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 T , we should integrate ODEs (40) and (41) twice with the initial conditions, i.e., 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 as Note that wave velocity c varies along the vessel approximately proportional to a −1∕2 , i.e. slower than the lumen radius. The variable varies fast only in a narrow and strongly tapering vessel. If we approximate c, and k = ∕c by their averaged values, we have Now ODE (48) admits an analytical solution in a number of cases. We consider here the case of a conic vessel segment: a(x) = a 1 + (a 2 − a 1 )(x∕L) , 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., where The transmission matrix for this vessel is A tapering vessel can be approximated by a sequence of truncated conic elements with values c, and k, approximated by their average value over every element. Let {x 0 = 0, x 1 , … , x N e −1 , x N e = L} be edge points of the cones. Here, N e is the number of conic elements.
The solution to the first initial value problem gives matrix elements T 11 and T 21 and the solution to the second initial value problem gives matrix elements T 12 and T 22 . This method is accurate and effective for the lowest harmonic components but looses its speed and accuracy for 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 x i . Using (52) we can write the transmission matrix T i of every element as Here, where a i is the lumen radius at x i . The transmission matrix T of total segment can be calculated through the product of transmission matrices T i for every element as: Accuracy of the piecewise conic approximation To study the accuracy of the piecewise conic approximation, we compute the reflection R( ) and transmission S( ) 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 L = 9.4 cm , the inlet and outlet diameters, respectively, are D 1 = 6.75 mm and D 2 = 3.5 mm . 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 h b = 2, 1 and 0.5 cm for the conic approximation. These lengths are adjusted to the vessel length by the rule: N e = ⌈L∕h b ⌉ , h = L∕N e , where ⌈x⌉ means the closest integer equal or greater than x. Therefore the element lengths are adjusted to values h ≈ 1.88, 0.94 and 0.49 cm.
The reflection R conic and transmission S conic coefficients computed in the conic approximation are compared with the same coefficients R numer , S numer computed by the numerical integration with the accuracy of 10 −10 . 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.

Boundary conditions
Consider a network containing N s vessels and vessel segments with one inlet segment s = 1 and a number of terminating segments. In every segment the axial coordinate x is referenced from its inlet. The flow parameters P ns (x) = P s (x; n ) and Q ns (x) = Q s (x; n ) in the sth segment can be calculated by the methods described above if the inlet flow parameters P ns (0) and Q ns (0) are known. Introduce the notations P ns = P ns (0) and Q ns = Q ns (0) , where s = 1, … , N s . Thus P ns ,Q ns are inlet amplitudes for every segment of the network. To compute the flow in the entire network, we have first to calculate 2N s unknowns [P n1 ,Q n1 , … ,P nN s ,Q nN s ] , for every harmonic component.
Hence, we have to derive 2N s equations with respect to these unknowns from the inlet, junction and outlet boundary conditions.
In Flores et al. (2016) the following 2N s unknowns are taken: [P n1 (0), P n1 (L 1 ), … , P nN s (0), P nN s (L N s )] . This is suitable if the vessels are approximated by uniform pipes without nonlinear correction. In stead of dealing with P ns (0) and Q ns (0) , we can actively use the transmission matrices T 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 p in (t) = p in (t + T) . We regard p in (t) as a full pressure taken from direct measurements (i.e. sum of forward and the backward waves). Then the boundary conditions is where P in n are calculated via the Fourier transform (35). If p in (t) 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 form where Ỹ n1 (0) and n1 (0) are inlet parameters of the first element in the first segment. Now eliminating the amplitude of the backward wave C b , we obtain an equation with respect to unknowns P n1 and Q n1 . Substituting P n1 and Q n1 , we have where Z n1 (0) is the characteristic impedance of very first element of the first segment.
Analogously, if the periodic flow rate waveform q in (t) = q in (t + T) is given at the inlet and it is the full waveform then the inlet boundary condition is where Q in n are calculated via the Fourier transform (35). If q in (t) is the only forward propagating waveform then for the nonzeroth harmonics, we can derive in the similar manner Matching conditions at junctions Consider the sth segment of length L s connected to  daughter segments at its outlet. Remind that P ns and Q ns are inlet flow parameters for the sth segment for the nth harmonic component. Now, the outlet (56) P n1 = P in n , n = 0, 1, … , N.
flow parameters can be written in the linear approximation through the transmission matrix T ns for the sth segment as Hereinafter T ns without argument means T ns = T ns (L s ) . 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 = 0 , i.e., If we denote the set of daughter segments associated with a parent segment s as = {d 1 , … , d  } , the continuity of the total pressure and the mass conservation at junction can, respectively, be written as Substituting relations (64) into conditions (61), we obtain  + 1 linear homogeneous equations as Matching 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 equations The number of equations for the junction always equals the nodal index of the junction:  +  . In the particular case of a junction with  = 2,  = 1 occurring among the cerebral arteries gives us three equations (61) P (1) ns (L s ) = T ns 11Pns + T ns 12Qns , and Q (1) ns (L s ) = T ns 21Pns + T ns 22Qns . (62) These equations are valid for zero harmonic too if the transmission matrix T 0p 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 n , the impedance load for a terminating segment s can be defined by single complex number Z out sn 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 R 1 and R 2 and a capacitor C, the impedance equals Z out ns = R 1 + R 2 ∕(1 + i n R 2 C). 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: Z out ns =Z ns = c s ∕(A s sn ) . 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 P f ∕Q f =Z∕(1 − i ) . Therefore the non-reflecting outlet impedance should be Here superscript 'out' denotes the value at the outlet of terminating segment s. Expressing P (1) ns (L s ) and Q (1) ns (L s ) through the transmission matrix T ns we obtain the linear equation for P (1) ns and Q (1) ns as Value Z out ns 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.
(70) Z out ns = P f ns (L s )∕Q f ns (L s ). Zeroth harmonic component The inlet/junction/outlet conditions should be considered separately for the zeroth harmonic component. The component P 0s 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 R out s can be set at all such segments as: where p out 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 = 3 . Some studies show that a value of = 3 is optimal for laminar flow, = 2.33 for turbulent flow and = 2.76 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)  At this point, we have N + 1 closed systems of 2N s linear equations with respect to P ns and Q ns for s = 1, … , N s and for all frequencies accounted: n ∶ n = 0, 1, … , N . Recollect that for a negative n we have the same solutions but complex conjugate. Then for monitored site x ∈ [0, L s ] of a selected vessel segment s, the waveform can be computed through the transmission matrix T s 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 P (1) n (0) and Q (1) n (0) 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, as Estimations 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 as After computation of n (2) (x, t) for the segment, we have to apply the Fourier transform (38) to relation (76) and integrate the inhomogeneous ODEs to give with zeroth initial conditions at the segment inlet: P (2) n (0) = 0 and Q (2) n (0) = 0 . 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 [P (1) n ,Q (1) n ] ) and the solution to inhomogeneous ODEs with zeroth initial condition (this addend just gives solution [P (2) n , Q (2) n ]). 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 as Subscript n is omitted here for brevity. To calculate entries of the T matrix, we use the conic approximation described in Sect. 4.4. To calculate the quadratures, we employ the Trapezoidal rule using cones' edge points {x 0 , … , x N e } as points of discretization. Now, the flow parameters at the outlet of the jth element are given as where Here T(x i ) = T i T i−1 ⋯ T 1 , where T i is given by formula (54). In the case of n = 0 , solutions to ODEs (80) and (81) are given as Again the quadratures can be computed by the Trapezoidal rule as

Nonlinear corrections to the junction matching conditions
With the account for nonlinear corrections, the outlet flow parameters take the following form (85) P ns (L s ) = T ns 11P ns + T ns 12Q ns + P (2) ns (L s ) where P (2) ns and Q (2) ns satisfy inhomogeneous ODEs (80)-(81) with zero initial conditions. Now, we can invoke the continuity of total pressure at a junction. Equation (85)  (1) ns + T ns

22Q
(1) ns . Now, substituting (87), we obtain  + 1 linear inhomogeneous equations as Analogous 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 (87) P nd (0) = P ns (L s 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 2N s algebraic equations with respect to P ns and Q ns is still linear, but most of the equations are inhomogeneous. After solving all N + 1 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 10 m , lumen radius of 1 cm and c 0 = 6.17 m/s , 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 q in (t) = 1 × exp{−(t − 0.05) 2 ∕0.01 2 } cm 3 /s. Non-reflection boundary condition is set at the outlet. The standard blood properties are used, i.e., = 1.04 g∕cm 3 and = 0.04 cm 2 /s. The velocity profile at inlet is taken Boileau et al. (2015) with = 9 that corresponds to = 11. 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 T = 4 s that is significantly greater than the time of pulse propagation along the pipe, t = 10∕6.17 = 1.62 s . The results of computation with the frequency independent = 11 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.
Using the approach developed here, this decay can be evaluated analytically. The wave number that accounts for viscosity is For the case under consideration = 2 × 0.14 = 0.88 s −1 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 gives i.e. all harmonic components with frequency much higher than are decaying at the same rate as Thus, 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 0.0713 m −1 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 = 11 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 = 4 (the value often taken in 1D arterial network modelling) is plotted using the dashed line in Fig. 4 (right). One can see that = 4 underestimates the pulse peak decay in the pipe with a 0 = 1 cm . This indicates that in a large artery like aorta = 11 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).

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 p , q . Equation (94) admits an analytical solution, for example, for the following linear wave where is an dimensionless amplitude factor: = p (1) max ∕ c 2 0 ; parameter determines the ratio of the constant and oscillating parts. Solutions to PDEs (94) satisfying boundary condition p (2) (0, t) =q (2) (0, t) = 0 are: where The 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 p (1) =q (1) = ( + (1 − ) cos( t − kx)), k = ∕c 0 (96) numerical computation is 0.06% for an element length of h = 2 cm. Solution (96)-(97) helps us to compare contribution of different nonlinear parameters in a waveform. The growing part of the solution contains the following factor The convection term, ( q 2 ∕A) x , mainly contributes to B + 1,2 and 2 = 8∕3 for the Poiseuille profile and the term A∕A 0 p x in Eq. (27) is unity in relation (98). The value of 1 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 1 . 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 c 0 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, T = 1 s , and the pulse duration, = T∕4 = 0.25 s. 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 q max = 0.1 cm 3 /s to ensure that the process is linear and = 4 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, q max . 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 p = p∕ c 2 0 and q = q∕(A 0 c 0 ) , which are proportional to the amplitude parameter, = q max ∕(A 0 c 0 ) . 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 q max = 10 cm 3 /s.
In Fig. 8 one can see the evaluation of nonlinear corrections p (2) and q (2) 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 p (2) ∕ 2 , q (2) ∕ 2 . 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 q max,ref = 10 −3 cm 3 /s, compute waveforms p ref and q ref and use it as a reference linear solution. Then, we calculate the waveforms p and q at a realistic amplitude q max and calculate the nonlinear correction by the formula, Note that the nonlinear distortion of the pulse is high in this network due to a long ( L = 1 m ) 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. (2) / 2

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 25 cm 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, D A , considered in Sazonov et al. (2017), i.e., 24 mm, 34 mm, 44 mm and 50 mm. All the pipes in the experimental setup, except one with the aneurysm have a constant diameter. Therefore, the calculation of the transmission matrix T 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: = 4∕3 , 1 = 1∕2 .
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 subjectspecific 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 a = √ A∕ is the lumen radius). Such approximations are mainly employed in 1D modelling, which can be written in the following forms: (101) c 2 = b 1 e − a + b 2 with b 1 = 188.7, b 2 = 21.19, = 9 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). 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 q in (t) is imposed at the inlet of the first segment as in Blanco et al. (2015); Boileau et al. (2015) with the heartbeat period of T = 1 s . All the tapering segments are approximated by truncated cones. We have set constant values = 4 and = 1 as in the computational model Mynard and Nithiarasu (2008). The constitutive relation is taken in the form of relation (9) (102) .10, 1 = 5.053, 2 = 0.1114 and thus the nonlinearity parameter 1 = 0.5 . The three-element Windkessel (lumped) model is employed as terminal boundary conditions. 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

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.  (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)  where u(r, x, t) is the axial velocity, v(r, x, t) is the radial velocity, p(r, x, t) is the pressure. Subscripts t, x, r denote partial derivative with respect to t, x, r, respectively. Boundary conditions at the pipe wall r = a(x, t) are: In a rigid pipe ( a ≡ const ), this boundary value problem admits exact solution for a harmonic perturbation: where p x is the amplitude of the pressure gradient; J 0 is the Bessel function of order 0 and is a complex viscous shear wave number. Any non-harmonic perturbation can be represented as a sum of harmonic perturbations. Note that solution (108) is not a wave like solution: the velocity is independent of x, i.e. the perturbation does not propagate along the pipe but appears in the same phase in the full pipe as soon as the pressure gradient is applied. If the pipe is not rigid, then a harmonic wave can propagate along the length, i.e., where k = ∕c is the wave number, c is the wave speed. Due to viscous losses, the wave is exponentially decaying; this decay is described by the imaginary part of k as where Imk < 0 . The pressure gradient is not constant anymore, but it varies harmonically in time and space, i.e., p x = ik p(r) exp {i t − ikx} where p = p − p 0 .
In an elastic pipe, the boundary value problem does not admit an exact solution; hence, we have to search for an approximate solution in the case of small wave amplitude and narrow pipe compared to the wavelength. Thus we introduce two small parameters: The first one is the Mach number-amplitude parameter and the second one is proportional to the axial to radial velocity (104) u t + uu x + vu r + −1 p x = u xx + r −1 (ru r ) r (105) v t + uv x + vv r + −1 p r = v xx + r −1 (rv) r r (106) u x + r −1 (rv) r = 0 (107) u(a) = 0, v(a) = a t . 110) [u(x, r, t), v(x, r, t) ratio. Here again subscript 0 refers to the basic, stress free state for a 0 , A 0 and to the wave velocity of small perturbation for c 0 . We now replace c and c 0 , respectively by c∕ and c 0 ∕ . This is equivalent to that of replacing k with k . We now represent the solution as a sum of Womersley solution and a small perturbation to it. We refer the variables of the Womersley solution with hat and the perturbation variables we refer with a tilde. Thus we search the solution in the following form: where functions ̂ and ̃ having dimension of velocity are introduced to simplify comparison of different terms in equations. In this solution, all the variables have the amplitude factor and the perturbation terms additionally have a factor . Representation of the pressure has a more complicated structure, necessary to balance the different order terms in the equations as it is shown below. Substituting solutions (112)-(114) into Eqs. (104)-(106) and keeping only the leading terms and omitting temporally the oscillation factor, we obtain where H.O.T. means the higher-order terms with respect to small parameters and . In Eq. (116), the leading terms are of the order of , whereas in (117) and (118) (120) v(r) = −ikp c 0 r 2 − J 1 ( r) J 0 ( a 0 ) .
Correction to the pressure p = c 0̃ can be calculated after substitution of solutions (120) into Eq. (117) and solving with respect to . This will give pressure variation in every cross section. The solution is cumbersome, and we do not present it here. Recollect that in the rigid wall case, the pressure is constant in every cross section.
Solutions (119) and (120) are still incomplete as the wavenumber k = ∕c is still unknown. Consider the mass conservation equation A t + q x = 0 . The flow rate in a cylindrical pipe can be calculated as q = 2 ∫ a 0 u(r)r d r . We write the mass conservation equation as: Differentiating the integrand and the upper limit we rewrite it as: In analogy with (115) we can represent the area as: and we rewrite (122) by substituting (112) and (115)   Using the Womersley solution, we can calculate for a harmonically oscillating flow with given frequency in cylindrical rigid and flexible pipes as (Fig. 13) The fastest way to compute the parameter is to use Padè approximation. The following approximation gives sufficient accuracy The maximal error is 2.5% at | | = 7.6.   (131)