Elastic jump propagation across a blood vessel junction

The theory of small-amplitude waves propagating across a blood vessel junction has been well established with linear analysis. In this study, we consider the propagation of large-amplitude, nonlinear waves (i.e. shocks and rarefactions) through a junction from a parent vessel into two (identical) daughter vessels using a combination of three approaches: numerical computations using a Godunov method with patching across the junction, analysis of a nonlinear Riemann problem in the neighbourhood of the junction and an analytical theory which extends the linear analysis to the following order in amplitude. A unified picture emerges: an abrupt (prescribed) increase in pressure at the inlet to the parent vessel generates a propagating shock wave along the parent vessel which interacts with the junction. For modest driving, this shock wave divides into propagating shock waves along the two daughter vessels and reflects a rarefaction wave back towards the inlet. However, for larger driving the reflected rarefaction wave becomes transcritical, generating an additional shock wave. Just beyond criticality this new shock wave has zero speed, pinned to the junction, but for further increases in driving this additional shock divides into two new propagating shock waves in the daughter vessels.


Introduction
The human circulatory system is formed of a bifurcating network of arteries and veins which conduct blood to and from the heart.Each heartbeat drives a rapid pulse wave through the system, which is easily detected in the larger arteries as a non-invasive clinical measure of cardiovascular health.Over the years there have been a wide variety of theoretical models to characterise pulse propagation though the cardiovascular system.Early models considered the propagation of small-amplitude pressure waves through a bifurcating network of elastic-walled tubes conveying fluid (see summaries in the monographs of Lighthill [1] and Pedley [2]).The inviscid flow in each constituent vessel is described through a characteristic impedance and the division of mass across each junction is characterised by a reflection coefficient; this simple description facilitates coupling into small networks [3].Such approaches have since been extended to larger networks to mimic the human cardiovascular system with physiological parameter values (e.g.[4,5]).This linearised analysis is typically sufficient to capture the dynamics because the flow remains strongly subcritical, where the maximal Mach number (the ratio of the maximal flow speed to the wavespeed) has been estimated as 0.2 − 0.25 in dogs, with humans believed to be similar [2].
However, abnormal physiological events can also trigger pressure wave propagation through the cardiovascular system.For example, incompetent aortic values can increase the volume of blood ejected by the heart, significantly enlarging the amplitude of the pressure pulse (sometimes known as a 'pistol shot' pulse) [6,7].Theoretical modelling of the propagation of these large-amplitude pressure waves is significantly more involved than the linearised analysis discussed above, where the models must now account for wave steepening and shock (or elastic jump) formation [8].The theory of shock waves propagating along a single elastic-walled blood vessel is similar to other classical hyperbolic systems, such as the dam break problem in shallow water theory [9,10].In particular, classical analysis using the tools of nonlinear waves [9,11] can predict the distance along the aorta from the heart where a shock will first form [2].
Large-amplitude pressure waves can also propagate as a result of external events.The present study is motivated by the aim to quantify the onset of retinal haemorrhage (i.e.bleeding of the retinal blood vessels at the back of the eye) following a traumatic brain injury, particularly in infants [12].One hypothesised mechanism is that retinal haemorrhage arises as a result of an abrupt increase in intracranial pressure (the so-called 'pressure-rise hypothesis') [13] which is then transmitted to the cerebrospinal fluid within the meningeal sheath of the optic nerve (the large cable of dendritic fibres connecting the eye to the brain) [14,15].Behind the eye, modelling predicts that such a pressure increase can drive a propagating shock wave along the nerve sheath, which is then significantly amplified by reflection at the blind end at the sclera and which can rupture bridging blood vessels and lead to optic nerve haemorrhage [16,17].However, the central retinal artery and vein traverse the optic nerve sheath as they enter the eye, and so are exposed to this acute increase in external pressure, which they then transmit to the retinal circulation.Modelling predicts that the amplitude and shape of the pressure wave transmitted into the retina is linked to the length of vessel confined within the optic nerve just before it enters the eye; for large enough perturbations this pressure wave can steepen to form a shock in the retinal circulation, where the tail of this incoming pressure wave is distributed over such long lengthscales that it can be approximated as spatially uniform [18].The clinical hypothesis attributes retinal haemorrhage to the propagation of this pressure wave, which will spread through the retinal circulation and drive abrupt expansion, and possible rupture, of more fragile vessels in distal parts of the network.Hence, an understanding of how largeamplitude waves propagate through blood vessel junctions is an essential component of any predictive model for retinal haemorrhage.As a precursor to developing large-scale computational models for the retinal circulation in response to a large-amplitude pressure perturbation, this paper focuses specifically on the propagation of a large-amplitude disturbance across a single blood vessel junction.
Within each individual blood vessel, the governing equations for incompressible blood flow are typically hyperbolic.Three-dimensional fluid-structure interaction models require careful numerical treatment to resolve the wave propagation (e.g.[19][20][21]).Such fluid-structure interaction approaches can be extended to consider three-dimensional flows through single blood vessel junctions (e.g.[22]) but extensions to model large portions of the cardiovascular system in this way are computationally demanding, particularly for networks composed of multiple generations (e.g.[23]).A more computationally tractable approach is to cross-sectionally average the three-dimensional governing equations and use a closure condition to derive a spatially one-dimensional representation of the flow.Even these reduced equations can present a significant numerical challenge to resolve the wave propagation, and are typically solved numerically using bespoke finite volume schemes (e.g.[24,25]).Such models have been used to investigate additional effects such as vessel tapering (e.g.[26]), wall viscoelasticity (e.g.[27,28]) and blood rheology (e.g.[29]).
These lower-order models for individual vessels can be extended to networks by coupling spatially one-dimensional flow within each vessels through boundary conditions imposed at single points which mimic the junction.These models impose conservation of (total or static) pressure and flow, either by directly imposing continuity across the junction using extrapolation along characteristics to inform the matching (e.g.[24,25,30,31]) or by introducing a separate junction compartment (e.g.[32][33][34]).
When the Reynolds number is large enough, these hyperbolic systems can exhibit wave steepening and thus the numerical schemes must be specifically developed to allow the propagation of shock waves.The seminal work on shock capture in blood vessels was conducted by Brook et al. [35], who developed a numerical scheme which exploits a local Riemann problem in each finite volume; this Riemann problem has subsequently been studied in detail using sophisticated constitutive laws for the vessel wall mechanics for both arteries [36] and veins [37].To extend this approach to circulatory networks the Riemann problem first needs to be adapted for use across a junction, known as the junction Riemann problem [38,39].In this case the junction boundary conditions are embedded at the point of discontinuity in initial conditions of the Riemann problem, similar to an abrupt discontinuity in vessel stiffness [36,40].The analysis is analogous to these classical Riemann problems while the flow remains subcritical, but extension to transcritical and supercritical flows is not straightforward [38].Subsequent analysis has posed iterative methods to deal with the supercritical regime [39], but what is currently missing is an understanding of how the flow can become transcritical, as may occur in physiological situations such as when a vein becomes highly collapsed [38].In the analysis below we extend the junction Riemann problem to include transcritical flow, which we show spans a wide range of the parameter space.
In this paper, we explore highly nonlinear wave propagation across a single blood vessel junction composed of a parent vessel bifurcating into two identical daughter vessels.For simplicity, the imposed pressure forcing at the inlet of the parent vessel mimics our prediction of the pressure wave transmitted into the retinal circulation following an abrupt rise in intracranial pressure [18].Within this setup, we revisit the junction Riemann problem and extend the existing analysis to include transcritical flow, demonstrating how resonance can produce additional shock waves in the system.These are similar in spirit to the resonances observed by Han et al. who considered the Riemann problem in a single blood vessel with an abrupt discontinuity in elastic stiffness [40][41][42] and analogous resonances found in shallow water theory [43,44].
The paper is arranged as follows.In Sec. 2 we summarise our numerical method for modelling the junction as well as a separate second-order analytical approximation, a continuation of the classical small-amplitude approach to the following order in perturbation amplitude.Predictions of this model are presented in Sec. 3, along with evidence of a limit point where the behaviour of the system changes fundamentally.In Sec. 4 we present our version of the junction Riemann problem, which we compare against our fully numerical model.Finally, in Sec. 5 we consider a vessel bifurcation with a jump in elastic stiffness, and in Sec.6 we consider the effect of underlying flow through the junction.We conclude with a discussion and an outline of further work in Sec. 7.

Model
Our model focuses on the propagation of a large-amplitude pressure wave through a single blood vessel junction which bifurcates symmetrically into two identical daughter vessels (see Fig. 1).In this study we ignore the complicated three-dimensional flows local to the junction induced by the geometry of the bifurcation, but instead focus on the global properties of the system induced by splitting the flow leaving the parent vessel and moving into the daughter vessels, as a pre-cursor to studying flow in large-scale cardiovascular networks.
Each constituent vessel is modelled as a long elastic-walled tube of length L * i (i = p, d), where subscripts i = p, d refer to quantities in the parent or daughter vessels, respectively.Note that the two daughter vessels are identical, so the single subscript d describes both simultaneously.Each vessel is parameterised by the spatial coordinate x * along the tube centre line.The junction is denoted as x * = 0, and so the parent vessel occupies −L * p ≤ x * ≤ 0 while the (identical) daughter vessels occupy 0 ≤ x * ≤ L * d .We systematically reduce the complexity of the three-dimensional flow along the constituent vessels, resulting in a model that has been used in many previous studies of blood flow in flexible-walled vessels e.g.[2,24,35,45], with a full derivation given in our previous study [18].The reduced flow equations are then spatially one-dimensional, where at time t the blood flow in each constituent vessel can be described by three quantities: fluid velocity along the axis of the vessel, u * i (x * , t * ); the tube cross-sectional area A * i (x * , t * ) and the local blood pressure p * i (x * , t * ), where i = p, d.Furthermore, each vessel has a (uniform) baseline cross-sectional area denoted A * i (i = p, d) and all vessels are subject to a uniform external pressure p * e .We characterise the elasticity of each vessel wall through an elastic stiffness parameter k * i (i = p, d), which can be estimated from the Young's modulus of the tissue and the wall thickness (see details in our previous study [18]).
The blood within each vessel is modelled as a homogeneous fluid of constant density ρ.Blood flow in the retinal circulation is viscous, but here the viscosity of the blood is neglected for two reasons.Firstly, this is the simplest version of the model which admits the formation of large-amplitude nonlinear waves (shocks and rarefactions).Secondly, the timescales of traumatic brain injury are typically so fast that viscous effects do not make a significant contribution (e.g.impact timescales for falls from altitude are on the order of milliseconds [46]).It would be straightforward to incorporate the effects of viscosity in the governing equations using an empirical term [18,35,47].
In this paper, the system is forced by a prescribed pressure profile at the inlet to the parent vessel with amplitude ∆p * , to mimic pressure wave transmission from the proximal end of the network (e.g. as a result of the transmission from acute fluid pressure increases in the optic nerve sheath [18]).For simplicity, at the outlet to the two daughter vessels a uniform fluid pressure is imposed by setting p * d (x * = L * , t * ) = p * 0 .The governing equations are non-dimensionalised by scaling relative to the properties of the parent vessel: we scale all lengths on the square root of the baseline cross-sectional area in the parent (A * p ) 1/2 , velocities on u * 0 = (k * p /ρ) 1/2 , time on t * 0 = L * p /u * 0 and pressures according to p * = k * p p + p * 0 .Note that dimensionless variables take the same symbol as their dimensional counterpart without the superscript * .This results in the following dimensionless groups for i = p, d denoting the aspect ratio of the parent tube, the dimensionless driving pressure, the dimensionless external pressure on the system, the dimensionless stiffness of each constituent vessel, the dimensionless length of each constituent vessel and the dimensionless baseline cross-sectional area of each vessel.Note that this choice of non-dimensionalisation gives k p = 1 and A p = 1, and in this study we restrict attention to the case where the parent and daughter vessels have equal base cross-sectional area (i.e.A d = 1).The governing equations of mass and momentum for each vessel, in the ideal limit as the Reynolds number approaches infinity, are given by where q i = u i A i is the fluid flux through each cross-section of the vessel (i = p, d).Note that the governing equations in the parent vessel hold on the domain −L p ≤ x ≤ 0, while in the daughter vessel they hold on 0 ≤ x ≤ L d .The system is closed using a constitutive law (or 'tube law') to mimic the wall elasticity in the form where Γ is a dimensionless function linking the cross-sectional area of each vessel to the local blood pressure [2,35].Specifically, we use a nonlinear tube law for the wall mechanics in the form where m, n ≥ 0 are parameters.This tube law is an algebraic function which resists both vessel expansion (for m > 0) and collapse (for n > 0).The parameters m and n can be chosen based on experimental measurements, although most data is obtained ex vivo.Here, two particular choices are used to illustrate our results: a nonlinear tube law where m = 10 and n = −3/2 (values proposed for the giraffe jugular vein [35] and often used for characterising human veins [18,31]) and a simple linear tube law where m = 1 and n = 0.This choice of tube law sets the nonlinear wave speed c i along each constituent vessel as This, in turn, is used to calculate the local Mach number, M i (x, t), at every spatial location along the vessel, a classical metric for describing transcritical flow in the form Our approach uses a one-dimensional model to describe the behaviour within each constituent vessel (Eqs.(2-4)), similar to many previous studies of cardiovascular flow.However, boundary conditions across the junction are required to couple these vessels together.In this case we impose conservation of mass and conservation of total pressure (e.g.[25]).Since the daughter vessels are identical it is assumed that mass of fluid leaving the parent is divided equally between them, which then gives two constraints where the parameter r i (i = p, d) represents the partition of mass across the junction so that r p = 1 and r d = 2; the choice of these constants can easily be generalised to capture junctions involving more vessels.We assume that the vessel cross-sectional area, pressure and flow velocity are spatially uniform in each constituent vessel before the pressure perturbation is applied at the upstream end.In most examples considered in this paper, in the initial configuration the vessels are all assumed to take their baseline cross-sectional area A i (x, 0) = Āi = 1 with no flow u i (x, 0) = 0 (i = p, d), so that p i (x, 0) = 0.These initial conditions automatically satisfy the junction boundary conditions (8)(9).
However, since arteries and veins have baseline flows in opposite directions, in Sec.6 below we consider initial conditions with underlying flow.To ensure this initial configuration is an equilibrium state it must obey the boundary conditions at the junction (8)(9).We prescribe the (spatially uniform) flow velocity in the daughter vessels as u d (x, 0), which implies that u p (x, 0) = (r d /r p )u d (x, 0) (given that in this study the baseline cross-sectional area of all constituent vessels are assumed equal).Note that this baseline flow must be driven by a prescribed upstream flux as it carries no pressure gradient in this inviscid model.
For simplicity, in this study we fix the domain lengths as L p = L d = 1 and set the external pressure as p e = 0.
In summary, this model is a simplified representation of blood flow across a single dichotomous junction, designed to capture the main features associated with large-amplitude wave propagation, and has been kept deliberately simple so we can focus on the dominant physical mechanisms and compare to the classical work involving small amplitude perturbations [1,2].
In what follows we characterise the system using a variety of measures, but focus particularly on the response of the system adjacent to the junction.We use a superscript (j) to denote a quantity evaluated at the junction when approaching from a given direction.For example, we denote A (j) We solve this theoretical model using three complementary approaches.Our primary approach is to use a numerical method based on the approach of Brook et al. [35], described in Sec.22.1.We further use an analytical approach which generalises the small-amplitude (linear) models described above, which is summarised in Sec.2.2 and Appendix A. Finally we consider our version of the junction Riemann problem in Sec. 4.

Finite-volume numerical method
Our primary method for solving this system involves a finite volume numerical method.For numerical convenience in these simulations the inlet pressure is abruptly ramped over a (short) timescale t ap (nondimensionalised on t * 0 ) from the baseline value (p p = 0) to the driving pressure p p = ∆p after which it remains constant, so that The system of mass and momentum equations (2-3) is solved numerically using a modified nonlinear Godunov method, where the spatial domains in the parent vessel (−L p ≤ x ≤ 0) and the daughter vessels (0 ≤ x ≤ L d ) are discretised onto a uniform mesh forming a line of finite volumes arranged in series.This method was first developed to describe unsteady flow in the giraffe jugular vein [35] and has subsequently been modified to analyse a line of vessels joined in series [18].At each time step, we calculate the cross-sectional area at the inlet of the parent vessel using the prescribed inlet pressure (10) via the tube law (4) and then compute the corresponding inlet flow velocity u p (x = −L p , t) by extrapolating the backward characteristics in the parent vessel assuming the flow is undisturbed far ahead.We proceed by solving a nonlinear Riemann problem for the cross-sectional area and flow velocity in each finite volume.
Note that the final mesh point in the parent vessel and the first mesh point in the daughter vessels are coincident at the junction (x = 0), though the corresponding flow variables can be different.We couple these variables by extrapolating the characteristics leaving each vessel and applying the junction boundary conditions (8)(9).For simplicity, zero pressure boundary conditions are imposed at the downstream end of the daughter vessels, but these have negligible effect on the predictions since perturbations do not reach the outlet within the simulation time.
We show below that this characteristic extrapolation technique is less robust at larger driving pressures where the rarefaction flow in the parent vessel becomes transcritical.In Sec. 4, we construct a nonlinear Riemann solver which explicitly accounts for the division of mass and any jump in material properties across the junction, showing that these predictions exhibit strong agreement with our numerical method.We will also compare our results at lower amplitudes to a modified linear theory described in the following section.

Modified linear perturbation theory
Our second approach to study this system is a modification of classical linear theory, which has been shown to approximate small amplitude waves propagating across a blood vessel junction [2].To improve this approximation for waves of slightly larger amplitude (the focus of this paper), in Appendix A this linear theory (for the linear tube law with m = 0 and n = 1) is extended to the following order in perturbation amplitude, calculating the (steady) amplitudes of the constituent waves.In response to an inlet perturbation to the cross-sectional area of amplitude ∆A (where ∆A ≪ 1) so that the cross-sectional area at the inlet of the parent vessel takes the form A in = A p (−L 1 , t) ≈ 1 + ∆A, this second-order theory predicts the amplitude of the reflected waves in the parent vessel adjacent to the junction as and the amplitude of the transmitted waves in the daughter vessels adjacent to the junction as ) 3 Baseline results In order to assess how a large-amplitude pressure wave divides across a symmetric junction, we first illustrate baseline examples using both the linear tube law (i.e.m = 1, n = 0, Fig. 2) and the nonlinear tube law (i.e.m = 10, n = 3/2, Fig. 3).In each case we plot several spatial profiles along both the parent and daughter vessels (Fig. 2a,3a) as well as a timetrace of the vessel cross-sectional area at selected points along the vessels (Fig. 2b,3b).An animation showing these dynamics is provided as online supplementary   material.In both cases the behaviour is qualitatively similar: the pressure increase at the inlet drives a rapid expansion of the entrance to the parent vessel (see timetraces in Fig. 2b, 3b blue solid line), triggering a pressure wave which propagates along the parent vessel towards the junction (Fig. 2a, Fig. 3a).This pressure wave steepens as it propagates, forming a nonlinear shock wave (or elastic jump): as the wave propagates the cross-sectional area abruptly increases at each spatial location (see time-traces in Fig. 2b, 3b, red solid line).Behind the elastic jump the vessel amplitude remains approximately constant, equal to that set by the driving pressure.Upon encountering the junction, part of the propagating pressure wave is transmitted into the two (identical) daughter vessels and part is reflected back along the parent towards the inlet (Fig. 2a, Fig. 3a red dashed line).This reflection is associated with a large increase in flow through the junction (Fig. 2aii, Fig. 3aii, red dashed line).The pressure wave transmitted into the two daughter vessels forms an elastic jump as it propagates, proceeding with almost constant amplitude towards the outlets (Fig. 2a, Fig. 3a, green dotted line).However, the reflected wave in the parent vessel takes the form of a rarefaction wave which propagates slowly back toward the system inlet.In this example the corresponding local Mach number of the flow is largest close to the junction but the flow always remains subcritical (Fig. 2aiii, Fig. 3aiii).In summary, in the examples the pressure perturbation at the inlet triggers a propagating elastic jump which is both reflected and transmitted at the junction, forming two new propagating elastic jumps in the daughter vessels and a reflected rarefaction wave in the parent vessel.The qualitative behaviour is independent of the choice of tube law.
In order to quantify the transmission and reflection as the pressure wave encounters the junction, Fig. 4 illustrates the amplitudes of the reflected and transmitted waves as a function of the driving pressure.As before, the behaviour is qualitatively similar for the linear tube law (Fig. 4i panels) and the nonlinear tube law (Fig. 4ii panels).As the driving pressure increases from zero, the amplitude of the elastic jump propagating along the parent increases accordingly, with a corresponding increase in the reflected and transmitted waves at the junction, evident in both the vessel cross-sectional areas (Fig. 4a) and the fluid flux through the junction (Fig. 4b), alongside an increase in the local Mach number around the junction (Fig. 4c).However as the input driving pressure increases further, the reflected cross-sectional area reaches a maximum and then starts to decline (Fig. 4a), despite a continued increase in the other amplitudes (Fig. 4a-c).In this regime the reflected Mach number in the parent vessel close to the junction increases rapidly for only small changes in the driving pressure, with the Mach number in the parent approaching unity for a finite value of the driving pressure (Fig. 4c); we denote this limit point as ∆p c .Beyond this critical driving pressure (i.e.∆p > ∆p c ) there is a noticeable change in behaviour.The local Mach number at the junction in the parent remains fixed at M (j) p = 1 at the boundary between subcritical and supercritical (Fig. 4c), the reflected rarefaction in the parent vessel remains pinned to the junction and reflected cross-sectional area in the parent increases almost linearly with the driving pressure (Fig. 4a).Interestingly, flow in the daughter vessels shows no indication of this limit point behaviour, where both the flow velocity and vessel cross-sectional area at the junction continue to increase with the driving pressure (Fig. 4).In summary, as the driving pressure increases the flow across the junction can eventually become transcritical, beyond which the reflected rarefaction wave in the parent vessel becomes pinned to the junction.
Numerical issues limit our ability to push this numerical method to larger driving pressures than those considered here.While the numerical method is stable for low (subcritical) driving pressures (i.e.∆p < ∆p c ), beyond the critical point ∆p c the numerical method exhibits time-dependent oscillations and/or extremely sharp gradients in the spatial profiles at the mesh points around the junction, consistent with previous observations for supercritical flows [32].
In order to validate these numerical observations, in Fig. 4(i panels) we compare our predictions with the linear tube law to those from the second-order analytical theory (Sec.22.2, Appendix A).As expected, we observe good agreement between the analytical and numerical approaches for low driving pressures.However, predictions of the two approaches diverge from one another well before the limit point.Our modified linear theory significantly improves the prediction of the reflected cross-sectional area at the junction for intermediate driving pressures (Fig. 4ai), at the expense of a slight reduction in agreement in other metrics (Fig. 4ci).However, this weakly nonlinear theory shows no evidence of transcriticality, and so to explain this transition we now proceed to derive a formal Riemann problem at the junction, analogous to the junction Riemann problems considered by others [38,39].This Riemann problem is derived in Sec. 4, and results are presented alongside a comparison to the predictions of the finite volume method.

Riemann problem at a symmetric junction
We construct the full Riemann problem in the neighbourhood of the junction, with the setup illustrated in Fig 5(a).In line with classical Riemann problems in flexible-walled vessels (e.g.[35,36]) we use the same nonlinear governing equations (2-4) but now prescribe a piece-wise uniform initial profile of both the cross-sectional area of the vessel and the fluid velocity (denoted by capital caligraphic letters), but with an imposed discontinuity in one or both of these quantities coincident with the junction (x = 0).Hence, we apply the initial condition Any disturbance will propagate outwards from the point of initial discontinuity, and at any later time t > 0 these initial conditions must still hold far enough away from the junction; these far-field boundary conditions replace the inlet and outlet boundary conditions considered in Sec. 3. Hence, upstream pressure and flow conditions are imposed through changes in the far-field conditions A p and U p .Note that the corresponding initial wavespeeds C i (i = p, d) can be calculated from Eq. ( 6).
As time progresses this initial discontinuity breaks up into left-and right-propagating waves with an intermediate region (often known as the 'star' region in shallow water problems) surrounding the junction where we must solve for the (spatially uniform) vessel cross-sectional area and fluid velocity.The division of fluid (and any jump in material properties) across the junction means that only two quantities can be conserved, which we choose to be the mass and total pressure (in line with boundary conditions (8,9)).Hence, this intermediate region must have a stationary wave at the junction (x = 0) to accommodate the jump in the other variables (this stationary wave is essentially a zero eigenvalue of the system when formulated as a matrix problem [36]).When expressed in terms of the model variables, we denote the piecewise uniform solution across this intermediate region as where x p (t) and x d (t) denote the edges of this intermediate region in the parent and daughter vessels, respectively; these edges are computed as part of the analysis below.The corresponding wavespeeds across the intermediate region, C (i = p, d) can be calculated from Eq. ( 6).Using the notation of the Riemann problem, our boundary conditions of mass and total pressure conservation (8-9) take the form where r p = 1 and r d = 2, as before.Linking this intermediate region to the two far-field regions (where the flow takes the initial values) are left-and right-propagating waves.In a classical Riemann problem these are either a rarefaction wave or a shock wave, leading to four possible responses when the bifurcation is symmetric (Sec.44.2), each a generalisation of the four standard responses in a single blood vessel [35,36].However, for flow in a flexible-walled vessel with an abrupt discontinuity in elastic stiffness, additional (resonant) flow structures are possible, as discussed by Han et al. (2015) [40,41].Two of these resonant cases prove relevant to the analysis of our symmetric junction system presented below, and emerge when the rarefaction wave in the parent vessel becomes transcritical, filling the entire space between the upstream far-field region and the stationary wave at the junction; this resonant interaction leads to at least one additional shock wave in the system.More details on these resonant cases are given in Secs.44.3 and 44.4 below.A more complete study of the parameter space, involving several more families of resonances, has recently been conducted [42].
We now discuss the the governing equations relevant to each case in turn.

Governing equations for shocks and rarefactions
The full governing equations (2-3) can be expressed in Riemann invariant form [2] where c i is the nonlinear wavespeed defined in Eq. ( 6) and R i± are the corresponding Riemann invariants of the flow (i = p, d).By definition, these Riemann invariants remain unchanged along the characteristic curves in the (x, t) plane governed by The discontinuous initial profile breaks up into constituent waves, which can either be a shock wave or a rarefaction, propagating into either the parent vessel or the daughter vessels.We now consider shock waves and rarefaction waves in turn.
For a shock wave propagating steadily along vessel i with velocity ṡi (separating the far-field region and the intermediate region), conservation of mass and conservation of momentum across the shock front [2,9] results in the constraints ṡi = ± A for i = p, d.The requirement for the shock to be propagating away from the junction means that for the parent vessel (i = p) the shock speed must be negative while for the daughter vessels (i = d) the shock speed must be positive.However, in both cases these two conditions combine to the same nonlinear constraint linking the far-field conditions to the flow in the intermediate region in the form (U for i = p, d.To be a physical solution this shock wave must satisfy the appropriate Lax entropy conditions [48], which for the parent vessel take the form and for the daughter vessels take the form In contrast, for a rarefacton wave linking the far-field region in the parent vessel to the junction region, the spatial profiles of the wave vary over a spatial interval.Any rarefaction adjacent to the farfield must be bounded by straight-line characteristics emanating from the origin of the x − t plane, and so a rarefaction in the parent vessel is confined to the interval (U p − C p ) t < x < U (j) p − C (j) p t, while any rarefaction in the daughter vessels is confined to the interval ( To fix the uniform values of the flow variables across the intermediate regions (on the junction side of the rarefaction) we follow along characteristics leaving the far-field region into the rarefaction toward the junction.The argument is standard [9,42], and it emerges that we can compute implicit equations for the spatial profiles across a rarefaction adjacent to the far-field in either the parent vessel or the daughter vessels.For example, for every x and t within the rarefaction fan in the parent vessel we can express the cross-sectional area profile A (r) p (x, t) as the solution of the implicit equation derived from the negative Riemann invariant Similarly, for every x and t within the rarefaction fan in the daughter vessel we can express the crosssectional area profile A (r) d (x, t) as the solution of the implicit equation derived from the positive Riemann invariant Simultaneously, we derive nonlinear constraints linking the far-field conditions to the intermediate variables by imposing continuity along the edge of the rarefaction meeting the intermediate region.These take the form These rarefaction waves cannot cross the stationary wave at the junction and must not outrun the far-field solution, and so for the solution to remain physical we impose the constraints The special case where the parent rarefaction becomes transcritical at the junction (i.e.U

Four standard responses consisting of single shocks and rarefactions
The standard responses in a classical Riemann problem consist of a single shock or rarefaction in each of the parent and daughter vessels.Allowing for every combination produces four possible cases.In these four standard examples we solve for four unknowns: based on the four far-field conditions (A p , U p A d and U d ).The two flow velocities at the junction can be expressed in terms of the corresponding cross-sectional areas using either Eq. ( 21) for a shock wave or Eq.(26,27) for a rarefaction wave.The junction cross-sectional areas are then calculated numerically using Newton's method by imposing the junction boundary conditions (16).Our strategy is to consider which solutions are physically meaningful by imposing Lax conditions (22,23) on shock waves and physicality conditions (28,29) on rarefaction waves.
Examples of the four standard cases are plotted in Fig. 5(b): rarefaction in the parent vessel and a shock in each daughter vessel (PRDS, Fig. 5bi); rarefaction in the parent vessel and a rarefaction in each daughter vessel (PRDR, Fig. 5bii) shock in the parent vessel and a shock in each daughter vessel (PSDS, Fig. 5biii); shock in the parent vessel and a rarefaction in each daughter vessel (PSDR, Fig. 5biv).The points in parameter space corresponding to these four examples are highlighted on Fig. 5(c).It emerges below that these four classical cases are insufficient to span the parameter space.Two further cases, where resonance between a rarefaction wave in the parent vessel and the stationary wave at the junction produces an additional shock wave are discussed below, where this new shock wave can be of zero speed pinned to the junction (Sec.4.3) or propagating into the daughter vessels (Sec.4.4).

Resonance with a zero speed shock: PR0DS
For the parameter values considered below, in PRDS the rarefaction wave in the parent vessel case can expand all the way to the junction and interact with the stationary wave, producing an additional zero speed shock wave at the junction; we denote this case PR0DS.This is analogous to Class C resonance reported by Han et al. (2015) for a vessel with an abrupt discontinuity in elastic stiffness [40].Here, the stationary wave at x = 0 splits into two stationary waves an asymptotically small distance ϵ ≪ 1 up and downstream of the junction, located at x = ±ϵ respectively.This division then introduces an additional resonant region around the junction.Following Han et al. (2015) [40], we expect that the stationary wave just inside the parent (at x = −ϵ) will be supercritical (i.e. the flows on either side will be supercritical) while the stationary wave just inside the daughter vessels (at x = ϵ) will be the subcritical (i.e. the flows on either side will be subcritical).
In the asymptotically small inner region between the two stationary waves (−ϵ ≤ x ≤ ϵ) we denote the flow variables with the superscript (s) .Across this region the material and environmental properties are spatially uniform, but take values intermediate to those in the parent and daughter vessels.In this study we allow for discontinuities in the vessel stiffness, where the intermediate value is denoted k s , and the mass distribution parameter across the junction, where the intermediate value is denoted r s .The vessel baseline cross-sectional area and external pressure can vary similarly, but here they are assumed equal in the daughter and parent vessels.We assume that these unknown properties between the stationary waves (i.e.those with a subscript s) can be scaled linearly between their values in the parent and daughter vessels using a single interpolation parameter s, such that s = 0 denotes the properties of the parent vessel and s = 1 denotes the properties of the daughter vessels; these intermediate values can be written In the resonant state, we immediately determine the flow variables on the parent side of the supercritical stationary wave by taking the outer limit (x → ϵ − ) of the rarefaction profiles across the parent (e.g.Eq. ( 24)), where the critical value of the parent cross-sectional area A (j) p satisfies the implicit equation with error O(ϵ), and it follows that in a similar manner to Eq. ( 21), these reduce to the single nonlinear constraint between the two pieces of the intermediate region in the daughter vessels (U To fully determine the PRDSS resonant profile we use Newton's method to solve the two conditions across the stationary wave (Eq.37) and the constraint across the resonant shock (Eq.39) coupled to the condition for the shock in the daughter vessel (Eq.21 with i = d) for the four unknowns A The pre-existing shock in the daughter vessels is subject to Lax conditions (Eq.23), whereas the new (resonant) shock is subject to Lax conditions based on the properties of the parent vessel (as this shock has originated from a resonance in the parent vessel) in the form A typical spatial profile for an example of PRDSS resonance is presented in Fig. 5(a,vi).The point in parameter space corresponding to this example is highlighted on Fig. 5(c).

The junction Riemann problem across the parameter space
To determine the solution of the Riemann problem at a given point in the parameter space spanned by the four far-field conditions (A p , U p , A d , U d ) our strategy is test for all six possible cases using a variety of initial conditions.We check each solution against their corresponding Lax and physicality conditions to determine if the solution is valid.Once a solution has been identified we use continuation methods to follow this solution branch until it terminates.In this way we construct large-scale sweeps of the parameter space and identify all feasible solutions.Across all cases tested we found no evidence of more than one solution at any point in the parameter space.However, our parameter searches were restricted to subcritical and transcritical flows, which are connected across the parameter space, and multi-valued solutions may become possible when the initial flow is strongly supercritical [36,49].
In order to assess how these six different families of solutions for the junction Riemann problem coexist, in Fig. 5(c) we plot an overview of the parameter space spanned by the far-field conditions in the parent vessel (i.e.A p and U p ) assuming that far-field conditions in the daughter vessels are held fixed (A d = 1, U d = 0); we use different coloured markers to indicate each family of solutions and compare the outcomes using the linear tube law (m = 1, n = 0 in Fig. 5c, left panel) to those from the more general nonlinear tube law (m = 10, n = 3/2 in Fig. 5c, right panel).In this first case we assume that all material and geometric properties of the constituent vessels are identical, and so the only discontinuity is due to the division of mass across the junction.Around the point of undisturbed flow (i.e.A p = A p = 1 and U p = 0) the system exhibits the four classical solutions in distinct (adjoining) wedges.Such a distribution of the classical solutions is also seen in the analysis of the Riemann problem in a single flexible-walled vessel [42].However, for larger values of the parent vessel far-field cross-sectional area, A p , the system adopts a PR0DS resonant solution.For the linear tube law the limiting value of A p for the onset of resonance varies reduces almost linearly with increasing far-field velocity in the parent U p (Fig. 5ci), while for the nonlinear tube law the limiting value of A p is almost entirely independent of U p (Fig. 5cii).Furthermore, for even larger values of the parent vessel far-field cross-sectional area, A p , the system transitions to PRDSS where in each case the limiting values of A p follow a similar trend to the onset of transcriticality.In summary, the division of mass across the junction from the parent vessel into the daughter vessels is sufficient to induce resonance for sufficiently large values of the far-field cross-sectional area in the parent vessel.
It emerges that if, instead, the far-field conditions in the daughter vessels are varied and those in the parent vessel are held fixed, the resonant solutions occur for small values of A d (not shown here, see [42]), but this regime is not the focus of the present study.
We note that the full numerical simulations in Sec. 3 exhibit spatially uniform profiles in the farfield (both ahead of and behind the propagating waves, see examples in Figs. 2, 3).This means that each simulation of the full numerical model effectively mimicks a Riemann problem at a single point in the parameter space spanned by the four far-field conditions (A p , U p , A d , U d ).For an initial shock wave propagating along the parent vessel driven by an inlet perturbation of amplitude ∆p, the far-field conditions in the parent must satisfy the constraints (19) separating the perturbed region behind and an undisturbed flow region ahead, so that Hence, as the driving pressure increases the solution of Eq. ( 41) traces a path through the parameter space spanned by the far-field conditions (A p , U p , A d , U d ).For example, for the bifurcation diagrams traced as a function of the driving pressure in Fig. 6, the corresponding path through the parameter space of the Riemann problem is shown as a thick black line in Figs.5(ci) and 5(cii) for the linear and nonlinear tube laws, respectively.In order to quantitatively compare the three approaches taken to study this system, Fig. 6 re-traces the bifurcation diagrams of Fig. 4 with the linear tube law (the nonlinear tube law is qualitatively similar) and plots the predictions of the full numerical simulations against the corresponding predictions of the Riemann problem and the analytical theory, showing the change in the vessel cross-sectional areas at the junction (Fig. 6a,c) and the flux through the junction (Fig. 6b,d).As the driving pressure increases from zero the Riemann problem exhibits a PRDS solution (consistent with the spatial profiles of the full numerical simulations shown in Figs. 2 and 3,) and the traces of the amplitude of the reflected and transmitted waves at the junction agree almost exactly with those from the full numerical simulations (Fig. 6); as noted previously, the analytical theories agree well for small driving, but eventually diverge as the driving amplitude increases.Furthermore, the solution of the Riemann problem becomes transcritical at ∆p = ∆p c ≈ 0.513 (for the linear tube law) (Fig. 6), almost exactly coincident with the predictions of the full numerical simulations, transitioning to a PR0DS solution.Beyond this limit point the predictions of the two approaches agree qualitatively, although there is an increasing discrepancy as ∆p increases (particularly in the flow through the daughter vessels, Fig. 6b,d) which we attribute to the numerical issues encountered in the full problem caused by an inability to resolve the zero speed shock at the junction predicted by the Riemann problem.For even larger values of the driving pressure the Riemann problem exhibits PRDSS, but the full numerical solutions do not converge in this regime; it is not clear if this lack of convergence is due to an inability to resolve the additional shock waves in the daughter vessels, or because the larger driving pressure imposes much larger spatial gradients which are difficult to resolve.In summary, this figure demonstrates that the predictions of the full numerical model and the Riemann problem agree very well in the subcritical regime, but also how the predictions gradually deviate from each other beyond the onset of transcriticality; the characteristic extrapolation technique used to match across the junction in the full simulations develops numerical issues, and is not capable of predicting the additional shock wave(s) predicted by the Riemann problem.
We now proceed to consider changes in vessel stiffness on either side of the junction (Sec.5) and an underlying base flow through the network (Sec.6), comparing each with the original numerical model presented in Sec. 3.

A discontinuity in vessel stiffness across the junction
Blood vessel stiffness changes gradually through the circulatory network, where vessels generally become more compliant with increasing distance from the heart.In order to determine how a change in vessel stiffness influences the predictions of the model, and in particular the onset of transcriticality, Fig. 7 examines the behaviour of the system with the linear tube law (the nonlinear tube law is qualitatively similar) when the daughter vessels are less stiff than the parent (Fig. 7i) and more stiff than the parent (Fig. 7ii panels).In each case we plot the parameter space of the Riemann problem spanned by the farfield properties in the parent vessel (Fig. 7a), isolate the path through the parameter space corresponding to the inlet perturbation from Eq. (41) (black line in Fig. 7a) and then compare the traces of the amplitude of the reflected wave in the parent and the transmitted wave in the daughters measured at the junction to those of the full numerical simulations (Fig. 7b).The parameter spaces are qualitatively similar for the two different values of daughter vessel stiffness, but the onset of resonance is pushed to larger values of  the parent cross-sectional area as the stiffness of the daughter vessels increases (Fig. 7).A similar trend is evident in the bifurcation diagrams, where the limit point for the onset of both PR0DS and PRDSS are pushed to larger values of the driving pressure as the vessel stiffness increases (Fig. 7b).Note that again we see increasing discrepancy between the full numerical simulations and the predictions of the Riemann problem once the driving pressure exceeds the limiting value (Fig. 7b).In summary, having the daughter vessels more compliant than the parent vessel promotes the onset of resonance, occurring for lower values of the driving pressure.

Base flow through the network
Cardiovascular networks exhibit an underlying base flow and so an arterial (venous) junction would experience net flow into the daughter vessels (parent vessel) over a cardiac cycle.In order to assess how this base flow influences the behaviour of the model with the nonlinear tube law (the linear tube law is qualitatively similar), Fig. 8 compares the predictions of the full numerical simulations to the predictions of the Riemann problem for both positive and negative base flow, again plotting bifurcation diagrams in terms of the driving pressure amplitude showing the amplitude of the reflected wave in the parent vessel at the junction (Fig. 8a) and the Mach number of the flow in the parent vessel at the junction (Fig. 8b).The behaviour in each case is qualitatively similar to the cases reported in Fig. 6.The critical driving pressure for the onset of resonance reduces as the base flow increases in speed towards the junction.Underlying flow towards the junction reduces the critical driving pressure, while underlying flow away from the junction increases the critical driving pressure.In summary, increasing the base flow through the junction promotes the onset of resonance through transition to transcritical flow.

Discussion
In this study we have considered the propagation of a large-amplitude pressure wave across a symmetric blood vessel junction (from a parent vessel into two identical daughter vessels) driven by an abrupt pressure perturbation at the inlet, as a precursor to studying the propagation of large-amplitude pressure waves in cardiovascular networks.We studied this system using three complementary approaches: nonlinear simulations using a generalisation of a specifically constructed finite-volume method with explicit shock capture [35]; a junction Riemann problem (e.g.[38]) and an analytical theory based on extending small-amplitude approaches (e.g.[1,2]) to include higher-order terms.
We have shown that for sufficiently low driving pressures the flow across the junction is subcritical, with all three approaches predicting that a single incoming shock wave is reflected as single rarefaction wave back along the parent vessel towards the inlet, and transmitted as shock waves in each of the two daughter vessels (Figs.2,3); predictions of the full numerical model and the junction Riemann problem agree almost perfectly in this regime, and the analytical theories asymptote to the others in the limit of zero driving (Fig. 4).However, as the driving pressure increases the spatial profile of the rarefaction wave in the parent vessel expands and eventually resonates with a stationary wave at the junction.This transition occurs as the flow in the parent vessel becomes transcritical i.e. the flow speed and the wave speed in the parent vessel become equal at the junction (Fig. 4c).Predictions of the critical driving pressure for the onset of resonance agree very well between the full numerical model and the junction Riemann problem, but this transition is missed entirely by the analytical theory (Fig. 6).Beyond this critical point the system enters a resonant state which exhibits at least one additional shock wave in the Riemann problem; predictions of the full numerical model and the junction Riemann problem agree qualitatively in this regime, though the full numerical simulations cannot capture these additional shock waves and exhibit some numerical issues, and for even larger driving the method eventually fails to converge entirely.An obvious item of future work would be to directly include our junction Riemann problem within the fully nonlinear finite volume code, replacing the characteristic extrapolation technique which is currently used to span the junction.Such a combined approach would then be capable of examining a much more general class of wave propagation problems, and would naturally extend to a large scale network.
Furthermore, we have shown that an abrupt decrease in vessel stiffness as the parent vessel divides across the junction (as expected in most areas of the human body) promotes the onset of resonance (Fig. 7).In addition, simple baseline flow driven through the junction (by a prescribed flow rate) was shown to promote (suppress) the onset of resonance when that baseline flow is directed toward (away from) the daughters as in arterial (venous) networks (Fig. 8).
We have shown how subcritical solutions of the junction Riemann problem can eventually become transcritical as the far-field conditions vary.However, we note that there may be other supercritical branches of solution entirely disconnected from the solutions identified here [39], facilitating regions of parameter space where the system has multiple co-existing solutions (e.g [36,49]).It should be noted that the range of far-field conditions considered here does span regimes where the initial conditions are supercritical.However, this supercriticality only plays a role in the solution for a very small wedge of the space (white region close to A p = 1, U p = 1 in Fig. 5c, left panel); connecting the supercritical solution in this wedge to the other solutions identified across the parameter space is an active area of future work.
The theoretical model adopted for this study is deliberately very simple, missing many of the features of human cardiovascular networks.In particular, the model employs only a very simplified description of the vessel wall mechanics, neglects vessel wall tapering and all viscous effects.The latter assumption is particularly limiting as the analysis not only misses the viscous losses associated with flow along a single vessel (which could be included using an empirical term in the governing equations [35,47]), it also neglects energy dissipated by the complicated three-dimensional flows generated in the neighbourhood of the junction (which can be computed e.g.[22]) and also the energy dissipated across the propagating shock fronts themselves as the vessel walls rapidly expand (see discussion in Pedley [2]).However, our model reflects the simplest version of the system which admits nonlinear shock waves, and it is an active area of future work to consider how the predictions are modified by viscous effects.
Funding: TAS and PSS acknowledge funding from EPSRC grant no.EP/P024270/1.ISO acknowledge funding from the Petroleum Technology Development Fund (PTDF).PSS acknowledges funding from EPSRC grant nos.EP/N014642/1, EP/S020950/1 and EP/S030875/1.( 9)) to obtain g (1)  p (z) ≡ g (1)  p H(z), g (1)  p = To improve this approximation we continue to the following order in perturbation amplitude.Following the same approach and now retaining terms of order ∆A 2 gives in each vessel (i = p, d) To the best of our knowledge this system cannot be solved analytically without further assumptions.To make analytical progress we assume that the local spatial derivatives are zero (which will be true away from the propagating wavefronts).In this case all forcing terms are zero and the second-order solution also follows a linear wave equation with wavespeed c i and so solutions can be written as (i = p, d) i (x + c i t) , q(2) i (x + c i t), in terms of wave profiles f i (z) and g i (z) in each vessel to be determined.Since there is no incoming perturbation at this order, all contributions are generated through interaction with the first-order flow at the junction and so f d = 0. Assuming the constituent waves are again square waves, applying the junction boundary conditions (Eqs.8, 9) the solution of the system takes the form g (2)  p (z) = g (2)  p H(z), g (2) d H(z), f These amplitudes are substituted back into the ansatz for the vessel cross-sectional area and flux (Eq.42) to obtain first-and second-order predictions for the amplitudes of the transmitted and reflected waves, as listed in Eq. (11), where k p has been set to 1 in accordance with our non-dimensionalisation.

Figure 1 :
Figure 1: Setup of the mathematical model with a parent blood vessel dividing into two identical daughter vessels: (left panel) a prescribed driving pressure applied at the proximal end of the parent vessel initiates wave propagation; (right panel) subsequent transmission of the pressure wave into the daughter vessels accompanied by partial reflection back along the parent vessel.

Figure 2 :
Figure 2: Spatial and temporal traces of the junction vessels using the full nonlinear model with a linear tube law m = 1, n = 0 in response to a pressure perturbation of amplitude ∆p = 0.5.(a) spatial profiles in the parent vessel (left) and daughter vessels (right) at t = 0.19 (blue solid line), t = 0.99 (red dotted line) and t = 1.49(green double dotted line), for (i) vessel cross-sectional area, (ii) flow velocity and (iii) Mach number; (b) the corresponding time-traces of the vessel cross-sectional area at three spatial locations along the parent vessel (the inlet, x = −L = −1, shown as a blue line, the midpoint, x = −0.5, shown as a red line, and adjacent to the junction, x = 0 − , shown as purple line) and two spatial locations in the daughter vessels (adjacent to the junction, x = 0 + , shown as orange line, and at the midpoint, x = 0.5, shown as green line).The spatial and temporal dynamics for this case are also shown in the accompanying video.

Figure 3 :
Figure 3: Spatial and temporal traces of the junction vessels using the full nonlinear model with a nonlinear tube law m = 10, n = 3/2 in response to a pressure perturbation of amplitude ∆p = 37.5.(a) spatial profiles in the parent vessel (left) and daughter vessels (right) at t = 0.046 (blue solid line), t = 0.093 (red dotted line) and t = 0.14 (green double dotted line), for (i) vessel cross-sectional area, (ii) flow velocity and (iii) Mach number; (b) the corresponding time-traces of the vessel cross-sectional area at three spatial locations along the parent vessel (the inlet, x = −L = −1, shown as a blue line, the midpoint, x = −0.5, shown as a red line, and adjacent to the junction, x = 0 − , shown as purple line) and two spatial locations in the daughter vessels (adjacent to the junction, x = 0 + , shown as orange line, and at the midpoint, x = 0.5, shown as green line).The spatial and temporal dynamics for this case are also shown in the accompanying video.

Figure 4 :
Figure 4: Predictions of the flow on either side of the junction as a function of dimensionless driving pressure using the (i) linear tube law (m = 1, n = 0) and (ii) the nonlinear tube law (m = 10, n = 3/2).Predictions of the fully nonlinear model show the amplitude of the transmitted wave (blue +, with points joined by a dashed line) and the reflected wave (black +, with points joined by a solid line) in terms of: (a) vessel cross-sectional area; (b) fluid flux and (c) Mach number.Predictions of the first-order (grey solid line) and second-order linear (grey dotted line) theories are shown for the linear tube law in (i).The dotted line in each panel indicates the point where the flow in the parent vessel becomes transcritical, i.e. ∆p = ∆p c .

Figure 5 :
Figure 5: Junction Riemann problem.(a) Notation for the junction Riemann problem (illustrated for a case with a rarefaction in the parent vessel and shock in the daughter vessels); (b) typical spatial profiles for the six different solutions to the junction Riemann problem (assuming a linear tube law, m = 1, n = 0) at the six points shown parameter space (c,i) below at time t = 0.25: (i) PRDS (triangle in (c,i)); (ii) PRDR (diamond in (c,i)); (iii) PSDS (circle in (c,i)); (iv) PSDR (cross in (c,i)); (v) PR0DS (square in (c,i)); (vi) PRDSS (inverted triangle in (c,i)); (c) Summary of the junction Riemann problem across the parameter space spanned by the far-field conditions in the parent vessel (A p , U p ) holding the far-field conditions in the daughter vessel fixed (A d = 1, U d = 0) for the (i) linear tube law (m = 1, n = 0) and (ii) nonlinear tube law (m = 10, n = −3/2), where the colour of the circle indicates the form of the response.In each panel in (b) the rarefaction waves are shown in blue and the shock waves are shown in red.

Figure 6 :
Figure 6: Predictions of the flow properties either side of the junction for increasing driving pressure using multiple methods.The linear tube law (m = 1, n = 0) is used here.(a) cross-sectional area in the parent vessel; (b) flow rate in the parent vessel; (c) cross-sectional area in the daughter vessels; (d) flow rate in the daughter vessel.In each panel the symbols are solutions from the full nonlinear computations, the solid lines are solutions of the junction Riemann problem where the line colour indicates the type of response and the grey solid (dashed) line is the prediction of the first-order (second-order) linear theory.

Figure 7 :
Figure 7: Effect of an abrupt discontinuity in vessel stiffness coincident with the junction.The linear tube law (m = 1, n = 0) is used here.(a) solutions of the junction Riemann problem across the parameter space spanned by the far-field conditions in the parent vessel (A p , U p ) holding the far-field conditions in the daughter vessel fixed (A d = 1, U d = 0) with (i) k d = 0.9 (ii) k d = 1.1;(b) tracing the properties of the flow on either side of the junction as a function of the dimensionless driving pressure comparing the full numerical simulations (symbols) with solutions for the junction Riemann problem for k d = 0.9 (+, solid lines) and k d = 1.1 (o and dotted lines) showing: (i) cross-sectional area in the parent vessel and (ii) cross-sectional area in the daughter vessels.The black lines in (a) trace the path through the parameter space of the junction Riemann problem which mimic the full nonlinear simulations, Eq. (41).

Figure 8 :
Figure 8: Effect of underlying base flow through the junction .The nonlinear tube law (m = 10, n = 3/2) is used here.The properties of the flow on the parent side of the junction are shown as a function of the dimensionless driving pressure comparing the full numerical simulations (symbols) with solutions for the junction Riemann problem for u p = 0 (+, solid lines), u p = 0.1 (o, dotted lines) and u p = −0.1 (×, dashed lines).(a) cross-sectional area in the parent vessel; (b) Mach number in the parent vessel .