TOWARDS THE DEVELOPMENT OF A FULLY COUPLED ARTERIAL-VENOUS 1 D MODEL : SUITABILITY OF USING A 1 D FINITE VOLUME METHOD WITH STAGGERED SPATIAL DISCRETIZATION

In this paper we outline the development of a 1D finite volume model to solve for blood flow through the arterial system. The model is based on a staggered spatial discretization which leads to a stable solution scheme. This scheme can accurately capture the various pressure wave reflections at locations with distinct discontinuities (in both area and material properties) as well as naturally treat branching vessels. We investigate the behaviour and performance of the solver for both a stented and branching vessel and finally for a full arterial network consisting of 55 arteries within the human vascular network.


INTRODUCTION
Within the cardiovascular system there are numerous examples of 'shunts'.In brief, a shunt, whether from a medical procedure or naturally occurring anomaly within the body, is an irregular passage of blood.In particular, we aim to investigate shunts between the arterial and venous system, for example a shunt medically inserted for the purposes of haemodialysis treatment, or anomalies such as the vein of Galen malformation.In both these instances, a shunt that allows for an abnormal amount of blood to pass from the arterial network back into the venous network often results in congestive heart failure.While several studies have been devoted to the investigation of flow and wave propagation through the arterial network, the venous system and the coupling between the two networks have received little attention.One dimensional transient models provide an attractive option for modelling the global flow phenomena through the arterial-venous network at computationally reasonable times.
Mathematical and numerical modelling of the cardiovascular system has received a lot of attention in recent years.Within this context, simplified models, including 1 dimensional models have demonstrated potential in providing useful information and insight at reasonable computational costs.Despite the reduced computational cost of 1D models, the cost can be excessive, especially considering a fully coupled arterial-venous network.This requires solving for a large vascular network constituting hundreds of branching vessels, solving for flow through the terminal coupling vessels (arterioles, capillaries and venules) as well as several stages of non-return valves within the venous network.With these computational demands in mind, it becomes necessary to have available an efficient numerical solution method.
To this end, we outline the development of a 1D finite volume solver discretized over a staggered grid.This is in contradiction to the published literature pertaining to the solution of blood flow through the arterial network.Most researchers have opted for solving the branching vascular tree through the use of discontinuous solution schemes, including the discontinuous Galerkin [6,7], locally conservative Galerkin [4] and finite volume Godunov schemes [1].While these solution schemes have been applied successfully to the branching vascular tree, they require the solution of an additional set of non-linear equations at each of the branching points.Consider Figure 1(a), a representation of a collocated discontinuous grid across a bifurcation.At the bifurcating location there are 6 unknowns that need to be solved.Therefore at each bifurcation, a non-linear set of equations relating mass flow rate, total pressure conservation and an additional 3 compatibility conditions based on the system characteristics need to be solved.
Using a staggered grid discretization (Figure 1(b)), there are as many equations as unknowns and can therefore be solved naturally with fewer degrees of freedom.The primary purpose of this paper is thus to analyse the suitability of using the FV method discretized over a staggered grid for solving a branching vascular network.We aim to investigate the stability and dissipation properties of the solution scheme.We will further demonstrate that the method is capable of solving for distinct discontinuities (in area and material properties) within an arterial vessel.We will lastly demonstrate the solution of bifurcations and solve the blood flow through a full arterial branching network.

GOVERNING EQUATIONS
The mass and momentum equations for a cylindrical vessel with compliant walls expressed in terms of area and velocity can be shown to be [5]: where A, p and u are cross-sectional luminal area of the vessel, fluid pressure and velocity respectively.K R is a function representing viscous losses.Assuming fully developed steady flow with a flat velocity profile the viscous losses can be approximated by where µ is the blood viscosity.
The mass (1) and momentum (2) equations however have 3 unknowns.To close the set of equations, an algebraic area-pressure relationship is used which is of the form: where p ext represents external pressure, A 0 is the initial unloaded area (artery area when p ext = 0) and β is a function describing material properties.In this paper we use a non-linear pressurearea relationship of the form [6] where β is defined as and h is the artery thickness, E is Young's modulus and σ is Poisson's ratio.

Staggered FV discretization
In this paper, equations ( 1) and ( 2) are solved using the FV method on a staggered grid (see Figure 2).This allows for as many equations as there are unknowns, and hence no additional compatibility equations are required.
Firstly, the mass and momentum conservation equations are modified to be in terms of the quantities to be conserved, i.e. volume flow rate Q = Au and total pressure p 0 = 1 2 ρu 2 + p.To this end, the mass (1) and momentum (2) equations become where Consider Figure 2(a).There are two sets of control volumes (CVs).Areas and pressures are associated with the one set of CVs, indicated by subscript i, and momentum conservation is associated with the other set of CVs indicated with subscript I.The midpoint of each of the groups is defined such that their quantities are at the face centres of the second group.
Firstly, let us consider the discretization of the momentum conservation equation (8).Under the assumptions of the FV method, we change equation ( 8) into its integral weak form over the control volumes depicted in Figure 2 and integrate over a time increment: where CV refers to the momentum conservation control volume.Using the divergence theorem: (10) becomes: Using backward differencing in time where superscript '0' refers to the previous time step t and no superscript refers to time step t + ∆t.Furthermore, using central differencing , (12) is rewritten as where I refers to the momentum conservation CV and i refers to mass conservation CV, with element length ∆x I , face area where integration in time is first order explicit if θ = 0, second order Crank-Nicholson for θ = 1 /2 and first order implicit for θ = 1.Assuming a uniform finite difference mesh (all ∆x I = ∆x i = ∆x ), and dividing through by A I and ∆t we obtain where the discretized total pressure is and

Discretization at branching vessels:
For a branching vessel we use a staggered grid shown in Figure 2(b), where area and pressure node values are located at the centre of the branching CV.With this grid discretization, both the momentum and mass conservation equations are naturally treated.
The mass and momentum discretized equations for the mesh in Figure 2(b) (for a uniform finite difference mesh where ∆x i = ∆x I = ∆x) at the branching point become: The discretized total pressure equations become In order to compute Q 2i at the bifurcation, we make the assumption that Based on this assumption, volume flow rate at the bifurcation can be approximated as It should be noted that when generating the mesh, the area and material properties at the bifurcation are interpolated to be consistent with this choice of Q interpolation.

Material and Area Discontinuities
Within the arterial and venous network there exists several discontinuities or sudden changes in both material properties (β ) and unloaded areas (A 0 ).These discontinuities can be due to cardiovascular diseases, for example plaque build up may lead to a decreased A 0 with an increase in β while the presence of an aneurysm increases A 0 while decreasing β .Furthermore, certain surgical interventions such as vascular stents or grafts may also introduce discontinuities in both A 0 and β .
For the staggered solution scheme, the discontinuities are defined in the centre of the momentum conservation CVs and then linearly interpolated to the face centres (see Figure 3).This is done during the initial mesh generation and has the effect of spreading the discontinuities across two CVs.This allows for the solution scheme to remain stable for large discontinuities and ensures that no spurious pressure waves are generated due to the presence of a discontinuity in either A 0 or β .

Boundary Conditions
Even though the discretized equations no longer explicitly require the system characteristics, they are still used when applying the system inlet and outlet boundary conditions.The system being solved constitutes a set of hyperbolic equations, and hence applying non-reflecting boundary conditions are non-trivial.By using the system characteristics, it is possible to apply non-reflecting inlet conditions (where the inlet flow remains unaffected by any backward running waves).Similarly, the outlet boundary conditions are based on a reflection coefficient [4], which allow for outlet boundary conditions varying from full reflection to no reflection.The derivation and implementation of the 1D system characteristics are well documented [4,5,6,7] and will therefore not be reproduced in this paper.

Solution algorithm
The discretized equations ( 16) through (18) constitute a set of coupled non-linear equations.In this paper we solve the discretized equations using the Gauss-Seidel algorithm [2].Gauss-Seidel is arguably one of the simpler methods for solving a set of non-linear equations.The decoupled equations ( 16) through (18) are solved iteratively per time step until convergence.The convergence criteria used is where ε = 1 × 10 −6 is used.N refers to the total number of elements and i + 1 and i to the current and previous solution iterations respectively.Gauss-Seidel typically requires a large number of convergence iterations, while each iteration is fairly cheap.The solution method does however perform well for the branching vascular problem, especially considering the small time steps necessary to maintain stability of the system.
It should however be noted, for large discontinuities within an artery segment (for example a stented vessel, see Section 3.1), Gauss-Seidel iterations may become unstable and fail to converge.To circumvent this instability it becomes necessary to reduce the time step size or make use of a relaxation method.One advantage of the Gauss-Seidel solution method, due to the decoupled nature of the solution of the 3 equations, is the simplicity with which to replace the pressure-area relationship, for example to a collapsible tube algebraic relationship for investigating flow through the venous network.
Unless otherwise stated, the problems investigated in this paper are based on the first order implicit (θ = 1) time integration.The use of implicit time integration allows the solution to remain stable for larger time steps while suffering from only minimal numerical dissipation.

Stented artery
The first test problem we consider is blood flow through a stented arterial vessel, first investigated by Quarteroni et al. [5].The vessel is depicted in Figure 4.The vessel has a length of 15cm, unit diameter with material properties β 0 = 451352dyne/cm 2 and a stented section of 5cm where β = κβ 0 and κ = 100.Blood density of ρ = 1g/cm 2 and viscosity of µ = 0.035 poise is used.The input pressure pulse is modelled as a half sine wave: .
The sine wave introduces a sharp discontinuous pressure gradient at the leading and trailing edges of the pressure pulse.
The pressure pulse through the stented artery is solved using 75 uniformly spaced elements with time steps ∆t = 1 × 10 −5 .The discontinuous material property jump is linearly interpolated over two elements (Figure 5(a)) as described in Section 2.1.3.
Figure 6(a) shows pressure-time plots at the 3 measurement locations (P, M, D) for both the healthy and stented artery.We directly compare to the digitised results obtained from Quarteroni et al. [5], based on the Taylor-Galerkin method, with 105 non-uniform elements and ∆t = 0.2 × 10 −6 .Quarteroni et al. modelled the stent discontinuity using a C1 continuous function (Figure 5(b)), where the function varies from β 0 to κβ 0 over a length of 1cm.
As a first level observation, the results in Figure 6(a) qualitatively exhibit the same trends.The pressure pulse propagates through the healthy artery undisturbed with visible dissipation due to the frictional losses.Furthermore, an increased proximal pressure (measurement point P) is observed due to partial backward reflection as the wave reaches the stented section.There is a further observable acceleration of the wave through the stiffened region.
There is however a marked phase shift between the results obtained using the staggered FV and Taylor-Galerkin schemes.The FV method further computes a lower proximal pressure peak for flow through the stented vessel.Both these observations can be quantified by the Taylor-Galerkin scheme of Quarteroni et al. resulting in lower pressure wave speeds.Comparing the wave speeds, a difference of approximately 15% is observed.From the system characteristics it can be shown that a finite amplitude pressure wave propagates through an artery at a velocity of u + c, where c is defined as [4] For physiologically real conditions c u.We can therefore reduce the wave speed by modifying the material properties such that β 0 = 0.7β 0 .The results for the reduced wave speed is shown in Figure 6(b) once again compared to the results obtained by Quarteroni et al.There is now a much better correspondence between the two set of results.The pertinent question thus remains which of the two schemes produces the correct wave speed.
To demonstrate that the wave speed resulting from the staggered FV scheme is indeed correct, consider Figure 7, depicting the pressure pulse propagation through a healthy artery for a modified material property β 0 = 317356.64dyne/cm 2.The modified β 0 relates to an unstressed wave speed c 0 of We analyse here the unstressed wave speed, since the peak wave speed of u + c is dependent on the solution of both primitive variables A and u. c 0 on the other hand is independent of the solution procedure, and can be seen as the analytical speed at which a wave travels through an initially unstressed artery.
Based on the modified material property (β 0 = 317356.64dyne/cm 2), the leading edge of the half sine pressure wave should reach the 3 measurement points P, M and D at 0.01, 0.02 and 0.03 seconds respectively (related to c 0 = 375 cm/s).The pressure vs. time plots in Figure 7 is computed using the same spatial and temporal discretization as was used for generating the results in Figure 6 (i.e.75 elements and ∆t = 1 × 10 −5 ) as well as for a refined discretization of 1200 elements and ∆t = 1 × 10 −6 .We further set µ = 0 such that any observable dissipation and dispersion is due to the numerical scheme only.
From Figure 7 it may be noted that both the coarse and refined discretized solutions result in the expected unstressed wave speed.The solution relating to the coarse discretization exhibits minor dissipation, which can more clearly be noted by the increasing dissipation of the leading edge of the sharp pressure gradient discontinuity imposed by the half sine wave.Despite the dissipative properties, the pulse velocities differ by less than 0.5% between the coarse and refined solutions.Furthermore, since there are no frictional losses (µ = 0) it may be noted that none of the wave amplitudes undergo a noticeable change as the waves travel through the length of the artery.We therefore conclude that the staggered, implicit in time, FV scheme is capable of accurately computing the wave speeds with only minor numerical dissipation.

Vessel branching
The arterial network can be viewed as the combination of a series of branching connections, where each branch is often not well matched.A non-matching bifurcation results in a partial backward reflection of a wave as the wave reaches these branching points.This gives rise to a large number of forward and backward running reflections through a branching arterial tree.It is therefore important that the solution scheme be capable of accurately capturing these reflecting waves.
Since the solution based on the staggered FV discretization does not explicitly rely on the system characteristics to treat the bifurcations, it is important to test whether the solutions obtained at the branching points remain consistent.To this end we analyse a bifurcating branch with an analytical linear reflection coefficient R b = 0. Sherwin et al [7] have demonstrated that the linear reflection coefficient at a bifurcation for a forward travelling wave can be written as: The linear reflection coefficient defined in equation ( 26) is strictly speaking not applicable for the non-linear model except for very small pressure pulses.To this end, we impose as input to the bifurcating branch a single half sine pressure pulse with unit amplitude and a period of 0.05 seconds where µ = 0 and ρ = 1.The pressure-time plots for the matching bifurcation is shown in Figure 9(a), where the measurement points are located at the centres of each of the respective branches.The pressure pulse from vessel A enters the two children vessels undisturbed with no spurious pressure oscillations.Furthermore, since there are no viscous losses, and both terminals B and C have equal areas and material properties, the full pressure pulse is transmitted into both children vessels.where the vessel properties are shown in Figure 8(b) relating to R b = 0.5.Furthermore, the terminal resistance of vessel B is set to 1 (i.e.fully backward reflecting at the outlet).Once again a pressure pulse with a period of 0.05 seconds is used to ensure that each of the reflecting waves are distinct.It can be noted that there is a 50% partial backward reflecting wave in vessel A as the pressure pulse reaches the bifurcation.The pressure wave in vessel B is fully reflected at the outlet, which upon travelling back to the bifurcation undergoes a negative reflection before partial transmission into parent vessel A and vessel C. The bifurcating problem is solved using 1200 uniform elements with step size ∆t = 1 × 10 −5 .A fine spatial discretization is required to resolve the sharp impulse like pressure wave introduced into the system.A pressure pulse period of 0.05 seconds is far shorter than typical pressure pulses produced by the heart but allows for the analyses of distinct pressure pulse propagations.

Arterial Network
The last problem we investigate is the pressure pulse propagation through 55 main arteries in the human vascular system.The arterial tree connectivity is shown in Figure 10 where the artery vessel properties are shown in Table 1.The vessel properties reproduced in Table 1 are adopted from Wang et al [8] and Sherwin et al. [6].The material properties were modified by Wang et al. from previous published data to ensure well matched bifurcations for forward travelling waves.This was done so as not to obscure any backward running reflections from the terminal vessels.The outflow conditions at the terminal vessels are treated as resistance elements, where the resistance for each terminal vessel is included in Table 1.
The cardiac input is modelled as a half sine pressure wave, with an amplitude of 15 × 10 3 dyne/cm 2 and period of 0.3 seconds, where the period of the full cardiac cycle is 1 second.Figure 11 depicts the input pressure pulse over 4 cardiac cycles, where the input wave is imposed as input to the ascending aorta (vessel 1).The density of blood is taken to be ρ = 1.021 g/cm 3  with viscosity of µ = 0.035 poise.
The branching network is solved with a total of 1850 non-uniform elements with implicit    The pressure pulse and velocity propagation measured at the entrance of the ascending aorta (1), left femoral artery (46) and the left anterior tibial artery (49) are depicted in Figure 12 for the first 4 cardiac cycles.The cardiac output compares well to published literature.It is however difficult to ascertain the physiological accuracy of the produced pressure and velocity outputs since the properties of the arterial vessels are still largely unknown.There are however a few key properties of the produced cardiac output that qualitatively compare well with in vivo observations.Firstly, there is a distinct formation of a diacrotic notch for the pressure pulse within the ascending aorta.The peak pressure pulses increase as blood propagates to the extremities of the arterial network, while the overall mean pressure decreases, and regions of large flow reversal are observed due to the interaction of backward travelling pressure waves.
As a final note, the terminal vessels were modelled using resistance elements.It has been demonstrated [4] that the capillary vessels exhibit properties of both resistance and capacitance.Future work will however focus on developing and coupling the venous network to the arterial network tree, where the flow through the arterioles, capillaries and venules will be modelled using an equivalent vessel formulation [3] and hence will no longer require the use of terminal vessel approximations.

Computational performance
To provide an indication of the computational time required by the solver, consider the results tabulated in Table 2.The reported computational times are for a single Intel Xeon 2.27GHz core using 3 different mesh sizes and 2 different time step sizes.
The average solution time per time step scales well for an increase in the number of system degrees of freedom, where the average computational time approximately doubles for a doubling in the number of degrees of freedom.Furthermore, despite the Courant number for the solution being greater than 1, it may be noted that the number of convergence iterations required by the Gauss-Seidel algorithm is relatively low, highlighting its suitability for this

CONCLUDING REMARKS AND FUTURE WORK
In this paper we have demonstrated that the finite volume method, discretized on a staggered grid is a viable solution methodology for solving the 1D vascular network.The solution scheme naturally treats bifurcations and vessel discontinuities.Future work will entail detailed investigations into alternative solution algorithms to solve the set of discretized equations as well as focus on the inclusion of the venous network.

Figure 1 .
Figure 1.Spatial discretization for a bifurcation using (a) discontinuous solution schemes, and (b) staggered solution schemes.

Figure 3 .
Figure 3. Illustration of linear interpolation of discontinuities (in A 0 and β ) across a cell from the centre of the momentum conservation CVs to the face centres.

Figure 6 .
Figure 6.Pressure pulse through the healthy and stented artery compared to digitised results of Quarteroni et al. [5], for (a) β 0 = 451352dyne/cm 2 as per problem definition, (b) modified β 0 = 0.7β 0 , to demonstrate that the wave speed obtained by Quarteroni et al. is fractionally too slow.

Figure 7 .Figure 8 .
Figure 7. Pressure wave propagation through a healthy artery with modified material properties β 0 = 317356.64dyne/cm 2for which the unstressed wave speed c 0 = 375cm/s.The leading edge of the pressure wave should therefore reach the measurement points P, M and D at 0.01, 0.02 and 0.03 seconds respectively.(a) Magnification of the leading edge of the pressure pulse, (b) full pressure pulse.

Figure 9 .
Figure 9. Pressure pulse propagation through a branching vessel for a (a) matching branching vessel (R b = 0), and (b) non matching branching vessel (R b =0.5).

Figure 11 .
Figure 11.Half sign pressure input pulse over 4 cardiac cycles.

Figure 12 .
Figure 12.Pressure and velocity measurements at the entrance to the ascending aorta, left femoral artery and left tibial artery.

Table 2 .
CPU time spent using a single 2.27GHz core for the 55 arterial network tree.