A Solution of the Junction Riemann Problem for 1D Hyperbolic Balance Laws in Networks including Supersonic Flow Conditions on Elastic Collapsible Tubes

The numerical modeling of one-dimensional (1D) domains joined by symmetric or asymmetric bifurcations or arbitrary junctions is still a challenge in the context of hyperbolic balance laws with application to flow in pipes, open channels or blood vessels, among others. The formulation of the Junction Riemann Problem (JRP) under subsonic conditions in 1D flow is clearly defined and solved by current methods, but they fail when sonic or supersonic conditions appear. Formulations coupling the 1D model for the vessels or pipes with other container-like formulations for junctions have been presented, requiring extra information such as assumed bulk mechanical properties and geometrical properties or the extension to more dimensions. To the best of our knowledge, in this work, the JRP is solved for the first time allowing solutions for all types of transitions and for any number of vessels, without requiring the definition of any extra information. The resulting JRP solver is theoretically well-founded, robust and simple, and returns the evolving state for the conserved variables in all vessels, allowing the use of any numerical method in the resolution of the inner cells used for the space-discretization of the vessels. The methodology of the proposed solver is presented in detail. The JRP solver is directly applicable if energy losses at the junctions are defined. Straightforward extension to other 1D hyperbolic flows can be performed.


Introduction
Numerical modeling of one-dimensional (1D) flow in networks of spatial domains joined by junctions offers a satisfactory compromise between the quality of the numerical predictions and the computational cost. There are a variety of 1D flow applications in the literature such as industrial piping networks, traffic flow, water flows in open channels or blood flow in the human circulatory system [1]. Flow physics at the bifurcation is complex in symmetric or asymmetric bifurcations and also in arbitrary junctions. One-dimensional modeling has proven to be an accurate tool in complex applications, such as computational hemodynamics in deformable arterial models under pulsatile flow, leading to competitive results if compared with three-dimensional (3D) simulations [2,3]. Flow in elastic and collapsible tubes is of great interest, and especially in the context of hemodynamics, physiology and medicine, and has been analyzed in many studies [4][5][6][7][8][9][10][11] in the last few decades.
However, the numerical modelling of 1D flow in networks is still challenging. In the 1D framework the junction is a singular point, where the numerical scheme cannot be directly applied and therefore internal boundary conditions must be prescribed [12], leading to the Junction Riemann Problem (JRP). A shortcoming of existing methods is their inability to deal with discontinuities or high subcritical, transcritical and supercritical flow through junctions, as in 1D channel networks, or to deal with transonic and supersonic flow conditions at junctions in physiological flows, such as in the venous system because of postural changes [1,13]. Existing methods are based on coupling approaches for energy or momentum conservation to the continuity equation and the characteristic equations [12,14,15] considering only subcritical or subsonic flow conditions. Those methods assume that the Riemann problem solutions involve only rarefaction waves.
The modeling of networks of vessels with highly compliant walls, such as veins, makes the problem even more challenging. When transmural pressure falls below a critical value in thin-walled elastic tubes, the vessel will collapse, and its cross-sectional area will diminish dramatically Flaherty et al. [16]. The fluid-structure interactions in veins may also include choking or limiting flow conditions, shock waves, and self-oscillation [8,17,18].
One alternative to overcoming the difficulties that the 1D approach has to face at the junctions, is to increase locally the number of spatial dimensions [1,19,20], avoiding in this way the computational cost of multidimensional modeling in the whole domain. In [1], channel networks were modeled with junctions where the two-dimensional (2D) shallow water equations were solved using the true geometry, whereas in channels, the usual 1D equations were solved. It is notorious how for transient gas flow modeling in [21], the flow inside the junction, defined as a 3D prismatic container, was modeled using the 2D Euler equations.
In the context of circulatory vessels, junction models can be defined using a confluence volume with assumed bulk mechanical properties [13,22]. The volume reference at zero transmural pressure depends on the grid size used to spatially discretize the 1D vessels and the stiffness of the confluence can be set using an average of the stiffness of the vessels connecting at the same junction. By combining the pressures at the junction and the neighboring cells it is possible to solve supersonic flow conditions [13] using a suitable numerical discretization.
In this work we will follow the approach based in the exact solution of the classical Riemann problem at the junction [9,10,12,23,24]. The exact solutions will be based on wave relations for shock and rarefactions (or compression and decompression waves) derived by a standard eigen-structure analysis of the hyperbolic system of equations for flow in collapsible tubes in [25] including transonic and supersonic conditions. The system will provide the wave relations that link the states provided by 1D domains sharing a junction, as in the case of the classic RP [26].
As pointed out in [24], by solving the JRP we want to obtain coupling conditions completely consistent with the underlying hyperbolic system of conservation laws, so the resulting JRP solver can be used in combination with any numerical scheme of arbitrary order of accuracy. The resulting states will be used to compute numerical fluxes needed by the explicit scheme to evolve the solution within the 1D domain. In this way, internal boundary conditions will be prescribed avoiding the use of container-like formulations for junctions or the use of multidimensional modeling.
The explicit scheme selected in this work to compute the updating fluxes in the cells within the 1D domain is the first order HLLS solver [27,28], based on the HLL numerical scheme developed by Harten et al. [29] and extended to deal with discontinuous source terms [27,30]. The first and third order in time and space HLLS was successfully applied in [28] to model flow in collapsible tubes with discontinuous mechanical properties including collapsed states including subsonic, transonic and supersonic flow under conditions of extreme vessel collapse and sonic blockage.
The numerical scheme in vessel networks must be completed with two types of boundary condition: external and internal (the JRP in our case). In both cases the time evolution of the system is governed by the state in the interior of the vessel and by waves which enter in the vessel from outside its boundary. Therefore, the definition of external boundary conditions in inflow or outflow sections is related to the equivalent problem for merging flows, splitting flows or any other possible combination in a confluence of vessels. In order to provide a suitable approach to the JRP, it is first necessary to explore boundary conditions at inflow or outflow sections, and then, extend the resulting conclusions to formulate the solution of the JRP.
The rest of this paper is structured as follows. In Sections 2.1 and 2.2 we present the underlying mathematical model and the numerical framework, based on the finite-volume method, respectively. The numerical scheme for vessels will be completed with external boundary conditions considering inflow-outflow sections of the network in Section 2.3. The solution of the JRP is presented in Section 2.4. First, in Section 3.1, we analyze the results given by the JRP solver by solving different RPs with discontinuous geometrical and mechanical properties and with exact solution. JRPs with three and four confluent vessels are presented in Sections 3.2 and 3.3 respectively. Final remarks are given in Section 4.

1D Mathematical Models in Arteries and Veins
Simplified 1D models in compliant vessels are used to represent the essential physical features of the propagation of pressure and flow waves. The resulting 1D differential equations for conservation of mass and momentum balance lead to a first-order and nonlinear hyperbolic system. The formulation is written as: with U = U(x, t) and where t is the time, x is the axial coordinate along the vessel, A is the local cross-sectional area, Q = Au is the volume flow rate with u the cross sectional average velocity in the axial direction, p(x, t) is the average internal pressure over the cross section expressing the response of the vessel wall deformation to pressure variations. Here, f is the friction force per unit length acting on the fluid in axial direction and the blood density is referred to as ρ.
A blunt velocity profile is assumed in this work by using κ = 1. Coordinate perpendicular to Earth surface, η, accounts for the gravitational forces due to the presence of gravity acceleration g.

Elastic Mechanical Properties of Vessels and Tube Laws
This system of equations is closed with a pressure-area relation of the form p(x, t) − p e (x, t) = p tr , where p e is the external pressure, and p tr is the elastic transmural pressure. The elastic transmural pressure depends on the vessel stiffness K o = K o (x), the reference area A o = A o (x) and the reference pressure p o = p o (x). Following [9], transmural pressure is assumed of the form with α = A/A o and σ the dimensionless transmural pressure difference. The vessel stiffness K o has different formulations for arteries, and exponents in σ are of the form m > 0 and n ∈ [−2, 0]. A o is the vessel cross-sectional area for which the transmural pressure p tr is p o . Transmural pressure for arteries is commonly defined setting m = 1/2 and n = 0, while for veins, m = 10 and n = −3/2 are used, as they are typical values for collapsible tubes and are derived from considerations made for the collapse of thin-walled elastic tubes. When a thin-walled tube collapses, there is a contact region of the internal vessel walls that divides the cross section into two tubes running in parallel.

Conservation-Law Form
An important feature of the model is that the system cannot be expressed in conservationlaw form when trying to extend the flux function F to involve the transmural pressure derivative in space [25]. However, if we can consider the vessel mechanical properties constant across a discontinuity, A o , c o and K o , we find that: with with c o the reference pulse wave velocity (PWV). In general, PWV can be expressed as: for a given vessel stiffness K o , or The system of equations can be now expressed using a conservation-law form with where the flux vector and vector of source terms that can be written in quasilinear form as follows: with where matrix J has two eigenvalues, λ 1 = u − c and λ 2 = u + c, and two real eigenvectors, e 1 = 1, λ 1 T and e 2 = 1, λ 2 T .

Numerical Method
All vessels of the domain, a 1D network, are discretized in computational cells of constant length size ∆x, leading to N k cells in each k vessel. In k vessel, its local numbering is i = 1, . . . , N k . We will call those cells involved in the coupling with internal (junctions) or external boundary conditions (located at the inflow-outflow sections of the network) terminal cells, where the local numbering is I k = 1 or I k = N k , with I k the index of the computational cell of vessel k linked to the boundary condition. Uniform initial conditions, U i = U n i will be defined at time t = 0 in all computational cells: The updating of the initial state in each cell, according to the Godunov first order method, is written as: with F ± i∓1/2 the numerical fluxes, with F = F(F , G χ ). Following the quasi-steady wavepropagation algorithm [31], the updating formula in (15) can be rewritten in terms of fluctuations, denoted here by δM. In this work, we define the following relation between fluxes and fluctuations as follows: with where in cases of exact equilibrium δM ± i+1/2 = 0. On the other hand, we will refer to the evolved states over time at the terminal interfaces in I k = 1 or I k = N k of the k-th vessel as U k = [A k , Q k ]. This state, generated by the imposition of external boundary conditions or by the solution of the JRP, will be used to compute the updating fluxes at its respective terminal interface: or in fluctuation form as In this work, we will extend the solution of the classical JRP in junctions for subsonic flow in order to deal with transonic and supersonic flow conditions.

Numerical Computation of Fluxes at the Inner Interfaces
When geometrical and mechanical properties along the vessel are not constant, the system of equations cannot be expressed using a conservation-law form. The transformation of the initial system of equations provides an equivalent system in conservation-law form, where an approximate flux function can be defined, allowing the application of the HLLS scheme. Numerical fluxes F ± i∓1/2 between each two adjacent cells i and i + 1 (inner cells of each vessel) will be computed using approximate solutions of the local RP posed for a time step t ∈ (0, ∆t), defined for the system in (1) as follows: With the independence of the approximate solver used, the Consistency Condition must be satisfied [32] when applied to the control volume represented in Figure 1. The control volume is limited by the maximum and minimum wave celerities of the system, −∆x/2 ≤ λ L ∆t and ∆x/2 ≥ λ R ∆t, with λ L ≥ λ R . Integrating (21) over the control volume the resulting condition is: where the δ symbol is used to express space difference, is written using an approximate JacobianÃ: involving the Roe average value: The source term is included in the Riemann solver as a singular source at the discontinuity point x = 0 [33]. Considering that source terms are not necessarily constant over time, the following time linearization is applied [33,34]: with A a suitable value of vessel area defined depending on the flow conditions [28] andf the integral of the friction force in the control volume. Using a discrete approach of the speed at the cell edge, c i+1/2 [28], the discrete source term G i+1/2 can be expressed as: with p d = p + gρη an equivalent pressure generated by the combination of the pressure and the gravitational force. The integral of the solution in the volume of interest in (22) is written as: with Matrix J is a local matrix that provides a set of two real eigenvalues λ 1 = u − c and λ 2 = u + c with λ 1 , and two eigenvectors: and can be presented as the discrete Jacobian matrix of a flux functionF e defined at each edge or interface between cells, that satisfies: Under this condition, the following approximate semi-discrete RP in conservation-law form is presented [28]: where the new flux functionF e is defined at each side of the RP as follows: allowing us to define the fluctuations as: involving a flux difference and a source term.

The HLLS Solver
The HLLS method can be formulated by appropriately imposing the Rankine-Hugoniot (RH) conditions over the slowest (λ L ) and fastest (λ R ) wave celerities of the RP in (32), and the RH relation for the steady contact wave, of speed λ 0 = 0, at x = 0, leading to the evolved fluxesF − e,i andF + e,i+1 at the left and the right sides of the RP [27,28]: In order to allow variations in the geometrical and mechanical properties of the vessel, the inter-cell fluxes F ± i+ 1 2 in the approximate Godunov method in (15) are given by a combination of functionsF e and F ϕ : leading to: In order to generate the intercell fluxes, it is necessary to compute the speeds λ R and λ L . Suitable approximations of λ R and λ L are discussed in [27,28]. Their analysis is not the focus of the present work.

External Boundary Conditions at the Inflow-Outflow Sections of the Network
The definition of external boundary conditions for flow at the inlet or outlet sections in single vessels is necessary prior the generalization to junctions for merging flows, splitting flows or other possible combinations in a confluence of vessels. In this section, we will focus on the resolution of the evolved state U k and its limits when defining boundary flow The following auxiliary function is defined depending on the position of the boundary cell [24]: We will distinguish between outflow sections, where flow at the boundary cell leaves the domain, (gu) n k > 0, and inflow sections, where flow initially enters at the boundary cell, (gu) n k < 0. We also define the following boundary Speed Index (SI): where when SI g < 1 subsonic flow conditions will be defined.

Non-Linear Waves: Shocks and Rarefactions
For any boundary cell, and depending on the sign of the jump between the evolved state and the initial state for α, dα k = α k − α n k , the solution will be a shock if dα k > 0, or a rarefaction if dα k < 0 [24,25]. To analyze the cases of interest and the possible limits, we define the celerities of the system as λ 1g = u k − g k c k and λ 2g = u k + g k c k . As the solutions can only be developed in the inner region of each vessel, we will only need to explore solutions where g k λ 1g ≤ 0, as shown in Figures 2 and 3, where wave patterns for rarefactions and shocks are plotted, respectively. Solutions will be analyzed considering whether they take place in an outflow or inflow section.  • In the case of a rarefaction in a k vessel, the solution is given by the characteristic field [24,25], and the jump condition across the rarefaction in invariant form becomes: with the following entropy condition: Considering that along the rarefaction waves the vessel area decreases, the equation in (42) can be written as or written as if expressed using variations in u. Both results suggest that velocity and flow variations must ensure: and, therefore, during the suction produced by a decompression wave, the vessel area decreases while flow and velocity increased/decreased for backward/forward waves, respectively. • If the solution in a k vessel is a shock of celerity S, that will travel in the −g k x direction, the Rankine-Hugoniot (RH) condition provides the solution. Solutions for right and left traveling shock waves can be defined using the following function [24,25]: with where the shock speed S k is given by: under the following entropy conditions: The RH condition applied to the mass conservation equation in combination with the direction of the advance of the shock or compression wave, leads to the following expression: and considering that f g > 0 in (48), we find that in boundary cells, Therefore, velocity and flow variations must ensure: and for compression waves, the vessel area increases while flow and velocity are decreased/increased for backward/forward waves, respectively.

Pulse Wave Velocity in Arteries and Veins
The parameters m and n in the dimensionless transmural pressure difference σ are different in arteries and veins and have an impact on the variation of the pulse wave velocity with the vessel deformation. Figure 4 shows the variation of the dimensionless PWV, c/c o , versus the vessel deformation α. While dimensionless PWV for arteries is a monotonically increasing function, dimensionless PWV for veins is a convex function with a minimum value of α min = (n/m) 2/(m−n) ≈ 0.72. This property has an interesting effect when observing the evolution of the solution for the SI in a rarefaction in veins, when the jump in the vessel deformation α includes α min . In this case, the maximum value of the SI appears when α = α min as shown in Figure 5.

Outflow and Inflow Boundary Conditions in the Subsonic Flow Regime
For any boundary cell, the solution in the subsonic flow regime is defined by the solution of the non-linear waves (conservation of the invariant in (42) for rarefactions and the RH condition in (48) for shocks) together with mass conservation at the boundary cells: with and only valid under subsonic flow conditions. The system of equations in (55) is non-linear and is solved here using Picard iterations between times t and t + ∆t. In terms of notation, iterations are denoted by superscript m. The solution at m iteration or the initial solution X m will be updated as follows: In order to find δX m , we first consider the Taylor expansions of the functions in E(X). If the error terms O i (δX 2 ) are assumed to be negligible, the system of equations can be written as follows: where J E (X) is an N x × N x Jacobian matrix defined over the function vector E(X). Now, all equations are expressed in linear form, and E(X m+1 ) = 0 is imposed and δX m is found to solve the equation The system of equations in (60) can be solved using the Gauss-Jordan elimination method, and δX m+1 is obtained. Once all variables are updated, we can check if the new numerical error in E m+1 is under a desired tolerance. Otherwise, the procedure is repeated.

Limitation of Boundary Conditions at Outflow Sections
Considering that we are interested in imposing boundary flow conditions Q bc (t) at boundary cells, it is useful to define the limits in the variations in the jump in flow (g (Q − Q n )) k , depending on the resulting wave and on the initial flow conditions, subsonic or supersonic. In this section, we first analyze the solutions when imposing outflow boundary conditions, with g k u n > 0.
• Case A1. Decompression waves under initial subsonic conditions, 0 < SI n g < 1 ( Figure 6). Suction is generated in the outflow section and a rarefaction or decompression wave appears according to (47). The value of the SI g increases in the g k x direction with |u | > |u n |. As the solution can be only developed in the inner region of the vessel, gλ 1g ≤ 0, the value of the speed index at the boundary, SI g , must be limited to conditions of sonic flow: Therefore, a critical limit in the velocity and flow, u c and Q c , appears when imposing the boundary condition: Then, if the limit in (61) is exceeded, the solution is given by the function F 1 (α) : and the solution for flow is given by Q bc = Q = Q c . • Case A2. Decompression waves under initial supersonic conditions, SI n g > 1 (Figure 7). For a supersonic initial state, both celerities λ 1g and λ 2g fall outside the domain of computation and the limit in (61) has already been exceeded. Further acceleration/deceleration would mean to force the development of the solution outside the vessel. Under these conditions no information can travel inside the domain and the initial state cannot change. We cannot specify any boundary conditions [35]. • Case B1. Compression waves under initial subsonic conditions, 0 < SI n g < 1 ( Figure 8). When the flow is decelerated/accelerated at the boundary in presence of backward/forward wave, the vessel area increases. The value of the SI g decreases in the g k x direction and sonic conditions cannot be achieved. • Case B2. Compression waves under initial supersonic conditions, SI n g > 1 (Figure 9). For a supersonic initial state, the wave λ n 1g falls outside the computational domain. When the flow is decelerated at the boundary, the shock must satisfy the entropy condition in (51), ensuring gλ 1g < 0. Therefore, sonic conditions cannot be achieved. Figure 6. Case A1. Flow acceleration/deceleration at an outflow boundary condition generated by a backward/forward rarefaction wave. Initial subsonic state. In this outflow section a rarefaction wave evolves from x max = λ n 1g ∆t to x min = λ 1g ∆t. Rarefaction is limited by λ 1g = 0 or x min = 0.   Figure 9. Case B2. Flow deceleration/acceleration generated by a backward/forward shock wave at an outflow boundary condition. Initial supersonic state. A shock wave moving inside the domain is generated if the imposed boundary condition is able to ensure gλ 1g ≤ 0 departing form gλ n 1g > 0.

Limitation of Boundary Conditions at Inflow Sections.
In this section, we analyze the limits and constrains that appear in presence of backward and forward compression or decompression waves for inflow sections with g k u n < 0. We define the following cases.
• Case C1. Compression waves under initial subsonic conditions, −1 < SI n g < 0 ( Figure 10). Under these conditions, the vessel area increases dα > 0 and the solution is a shock or compression wave moving in the −g k x direction, where λ 1g,k > λ n 1g,k . • Case C2. Compression waves under initial supersonic conditions, SI n g < −1 ( Figure 11). In the supersonic range, both λ n 1g and λ n 2g fall inside the domain. Apart from the prescribed flow, another condition, such as vessel area or total energy must be imposed. Initial conditions and boundary conditions may satisfy the relation in (51) and entropy conditions in (48). • Case D1. Decompression waves under initial subsonic conditions, −1 < SI n g < 0 ( Figure 12). In the subsonic range, we find that |u | < |u n |. The limiting value of Q bc is given by g k λ 1g = 0, or g k u = c , and represents a condition of reversal flow. When this limit is exceeded, sonic conditions must be imposed, and the solution is given by function in (63). • Case D2. Decompression waves under initial supersonic conditions, SI n g < −1 ( Figure 13). Again, the limiting state, a condition of reversal flow, is provided by g k λ 1g = 0 and, if exceeded, the solution is given by function in (63). Figure 10. Case C1. Flow deceleration/acceleration, g k dQ bc < 0, at an inflow boundary condition under initial subsonic conditions. In this inflow section a backward/forward shock wave moving inside the domain is generated. Figure 11. Case C2. Flow deceleration/acceleration, g k dQ bc < 0, at an inflow boundary condition under initial supersonic conditions. In this inflow section a backward/forward shock wave moving inside the domain is generated. Both λ n 1g and λ n 2g fall inside of the domain. Two boundary conditions are required. Figure 12. Case D1. Flow acceleration/deceleration, g k dQ bc > 0, at an inflow boundary condition under initial subsonic conditions. In this inflow section a backward/forward rarefaction wave appears. The solution evolves from x max = λ n 1g ∆t to x min = λ 1g ∆t. When λ 1g = 0 the flow velocity direction is opposite to the initial one. Figure 13. Case D2. Flow acceleration/deceleration, g k dQ bc > 0, at an inflow boundary condition under initial supersonic conditions. In this inflow section a backward/forward rarefaction wave can defined inside the domain as λ n 1g falls within the domain. When λ 1g = 0 the flow velocity direction is opposite to the initial one.

Outflow and Inflow Boundary Conditions including Subsonic, Sonic and Supersonic Flow
The previous cases in Sections 2.3.4 and 2.3.5 can be reorganized depending on the initial flow regime. For subsonic flow, SI n g < 1, we find that: • If 0 < SI n g < 1, at outflow boundary conditions: 1.
The acceleration/deceleration of the flow generated by a backward/forward wave can be imposed provided that the limiting value of flow in (62) is not exceeded, leading to a decompression wave (case A1).

2.
The deceleration/acceleration of the flow generated by a backward/forward wave can be imposed leading to a compression wave (case B1), and no sonic limitation appears.
Flow deceleration/acceleration, generated by a backward/forward wave can be imposed leading to a compression wave (case C1).

2.
Flow acceleration/deceleration, generated by a backward/forward wave leads to a decompression wave. Flow at the boundary must by limited by the function in (63) (case D1).
For supersonic flow: • If SI n g > 1, at outflow boundary conditions: 1.
The acceleration/deceleration of the flow cannot be imposed, as no information can arrive to the inner domain of the vessel (case A2), The deceleration/acceleration of the flow, generated by a backward/forward wave can be imposed if leading to a compression wave with a suitable value of celerity (case B2), −g k S > 0.
The deceleration/acceleration of the flow generated by a backward/forward wave can be imposed, leading to a compression wave, provided that both flow and area are imposed (case C2), The deceleration of the flow can be imposed, leading to a decompression wave (case D2) limited by the function in (63).
Therefore, we need to provide solutions considering that: • The target functions in system (55) must be redefined to ensure limitation of sonic flow in cases A1, D1 y D2. • Unfeasible acceleration/deceleration, in supersonic outflow sections in case A2 must be avoided, the target functions in system (55) must be redefined to include this case. • In case C2, involving a supersonic inflow section, both area and flow rate must be imposed. In this case, the system of equations in (55) is useless.
In order to limit the possible transitions, the system of equations in (55) is modified. Now, error function in (55) is written as the sum of two functions: The first equation in (65), accounting for the solution of non-linear waves, involves the following velocity and δα = (α − α n ) k . T w is zero when sonic conditions are reached during the iterative search of a rarefaction solution, resulting in R w = 1. In this case, u is no longer an unknown and the rarefaction waves remain limited by λ 1g = 0. The coefficient D w retains the initial values for velocity and area for those outflow sections where supersonic conditions are initially defined (case A2) when necessary. D w is defined as: avoiding further acceleration of the flow in supersonic conditions or unphysical shock wave directions.
Mass conservation in the system of equations in (65) is defined using the following auxiliary variable leading to the initial value of discharge when D w = 0, to conditions of sonic flow when R w = 1 or to the objective value Q m bc,k if subsonic conditions are retained, T w D w = 1. The auxiliary flow value Q m bc approximates Q bc (t) and ensures a smooth convergence to the solution, will be defined below. Matrix In this way, when sonic conditions are reached and u is no longer an unknown, the numerical method presented here avoids the updating of u , setting δu = 0.

Initial Values and Relaxation Parameters
In the first iteration, m = 0, the solution will be defined assuming that: • if g(Q bc − Q n ) ≥ 0, the solution will be a rarefaction, • if g(Q bc − Q n ) < 0, the solution will be a shock.
In order to avoid an excessive variation of α between two steps, we reduce the step size δX m by using a damping factor γ α ≤ 1: with δα lim a constant value, defined to limit the variation of α. Next, we also limit the variation in δX m by imposing a constraint in the maximum variation of the velocity u, δu lim , by using the damping factor γ u ≤ 1.
In this case, δu lim is not a fixed value as δα lim . The variation in δu lim is computed assuming that the increase in SI g from iteration m to m + 1, δSI m g = SI m+1 g − SI m g , is small and small enough to control the solution when the solution is in the vicinity of sonic flow conditions. This is ensured by limiting the increment δSI m g by a prefixed constant value, δSI lim , which can be also reduced to avoid sonic conditions when SI g = 1, as follows: With this value of δSI m g , the maximum variation of velocity is approximated by: as small variations in α have been enforced. The final solution will also be achieved not directly imposing the value of Q bc in the mass conservation equation, but changing its value smoothly using the auxiliary variables Q m bc and δQ m bc , with: untill the final solution for Q bc is achieved, with: limiting the value of Q bc in the mass conservation equation when sonic conditions are reached. The details of the invariant integration are given in Appendix A.

Solution of the JRP
The Riemann problem at vertex P or JRP, for N P infinitely long converging vessels sharing vertex P reads: for k = 1, ..., N P . To couple the N P vessels, we have to compute the unknown cross-sectional area A k and velocity u k for each vessel converging at node P, which means that we have 2N P unknowns. Quantities A k and u k will be used to compute fluxes at the boundary interface of the k-th vessel. For the N P boundary cells meeting at the junction, the following mass conservation equation is written: where Q k = u k A k . The solution of the JRP under subsonic flow conditions is next revisited and later extended to general flow conditions.

Solution of the JRP under Subsonic Flow Conditions
The solution of the JRP under subsonic flow conditions can be defined by a system of 2N P equations of the form E(X) = 0. The first N P equations are relations for nonlinear waves (conservation of the Riemann Invariants in (42) or shock waves in (48)). The remaining N P equations, mass conservation in (80) and the N P − 1 equations involving conservation of total pressure, p T , are related to the stationary contact discontinuity generated by variable mechanical and geometrical properties [24]: (57), where vessel k = 1 is used to define the reference value in the energy conservation equation.
The system of equations in (81) can be linearized and solved using Picard iterations between times t and t + ∆t, using a Jacobian matrix J m E defined over the function vector E, where ffiX m is obtained. In this case, matrix J E is given by: where and with ∂σ ∂α

Solution of the JRP under Subsonic and Supersonic Flow Conditions
The system of equations for subsonic flow conditions is modified here to deal with both subsonic and supersonic flow conditions. After each iteration, the variables of the system in X and a set of dynamic parameters for each vessel will be updated. The new system of equations is split into two terms: is defined to impose conservation of the initial value conditions in a k vessel depending on the dynamic parameter D w,k . Parameter D w,k is defined as follows: It is worth noting that, for outflow sections under supersonic conditions, it becomes zero if the iterative solution leads to an unphysical solution. Unphysical solutions are associated with flow acceleration, to a reduction of the vessel area or to a shock wave traveling in the wrong direction (analogous to case A2 in the external boundary analysis). In this case, both the area and the velocity of the vessel are no longer an unknown of the system, and U k = U n k must be enforced. On the other hand, when D w,k = 1, the transmission of a backward wave generated downstream of a convergence point is guaranteed.
When coefficient D w,k = 1, the error function E k for vessels k = 1, . . . , N p , reduces to E a k accounting for relations across non-linear waves, involving the auxiliary variable u • in (66). If sonic flow conditions are reached during the iterative process for vessel k, the evolved solution for the velocity is g k c , and this velocity is no longer an unknown of the problem. When the parameter D w,k = 0, the error function E k = E b k , expressing that the initial condition for the velocity must be kept over time, u k = u n k . The definition of the velocity in (66) enables a generic formulation for the flow in a vessel as follows: and error function E N P +1 for mass conservation is split in E a N P +1 and E b N P +1 depending on the value of D w,k .
Error functions E N P +1+k for k = 2, . . . , N p are modified to consider that the evolved state of a vessel may not participate in the energy conservation equation. The reference value in the energy conservation equation is: with Coefficient L w,k is always zero except for in the first eligible vessel where energy conservation can be defined and this k vessel will be used as the reference value in the energy conservation equation. Ineligible vessels do not participate in the energy conservation equation when T w,k = 0 (flow limitation) or D w = 0 (conservation of the initial condition). Coefficient Q w,k is defined as: and may become equal to 0 or 1. In the case that L w,k = 1, Q w,k becomes zero as this k vessel will be used as the reference value in the energy conservation equation. Coefficient Q w,k also becomes zero when the flow reaches sonic conditions in a rarefaction or when the flow is supersonic and cannot be further accelerated with independence of the gradient of energy across the junction. In both cases, these ineligible vessels do not participate in energy conservation equation, as their solution becomes independent of the energy level of the other vessels. When for k > 1, D w,k = 0, the energy conservation equation in E a N P +1+k is omitted and preservation of initial conditions for vessel area, , with k the first eligible vessel with L w,k = 1. This k vessel does not require further equations as it has already been used to provide the energy reference level.
System of equations in (87) is indirectly solved using the equivalent linear system in (60), where Jacobian matrix J E (X) is defined as: and matrices J E, 1 and J E,2 involve relations for non-linear waves and matrices J E, 3 and J E,4 are related to the mass conservation and energy conservation equations. The Jacobian matrix J E (X) must be carefully constructed to guarantee that the system of equations provides a suitable solution for ffiX m when ineligible vessels appear, as the number of unknowns may change along the iterative process. When presenting the new system of equations in (87), instead of dynamically changing the size of the system of equations, we have enforced the necessary conditions to retain the original dimensions of the problem with the help of the dynamic coefficients. The same idea is applied to the Jacobian matrix J E (X).
In general, the number of unknowns and original equations is reduced to . Therefore, the Jacobian will have to deal with two situations, D w,k = 0 and R w,k = 1. When D w,k = 0 for vessel k, the linearized system will ensure that U k = U n k by returning δu k = 0 and δα k = 0. When R w,k = 1 for vessel k, the linearized system will return δu k = 0, as the solution will be provided by α .
Matrices J E,1 and J E,2 are given by: with When sonic conditions are reached in a k vessel, R w,k = 1, coefficient e k,k = 0, and δu k is not any longer an unknown of the system. If D w,k = 0, e k,k = 1 and e k,N P +k = 0, enforcing 3 is of the form: where for those ineligible vessels with D w,k T w,k = 0, coefficients e N P +1,k to e 2N P ,k become zero. Matrix J b E,3 is given by: and is defined to return δu k = 0 for those vessels where sonic conditions are reached. If D w,k = 0, the k column in matrix J b F,3 becomes null and vessel k does not participate in the energy error functions. Matrix with and, when Q w,k = 0, vessel k does not participate in the energy error functions. The same applies to matrix J b E,4 defined as is defined to return δα k = (α n − α ) k for those outflow sections where supersonic conditions are initially imposed and the solution will not change over time.
The system of equations is solved using the Gauss-Jordan elimination method, and δX m+1 is obtained. The updating step δX m is corrected by using a damping factor selected as the minimum among the damping factors γ α,k in (73), defined for each k vessel of the junction: Next, the updating step δX m is corrected again by using another damping factor, this time selected as the minimum among the damping factors γ u,k in (74), where sonic flow conditions are controlled by computing γ u,k as in (76).
Once the new value of X is computed, the values of R w,k , T w,k and D w,k are updated, and the velocity and vessel deformation are defined as For each k vessel where initial supersonic conditions are imposed, SI n g ≥ 1, only flow variations are allowed when D w,k = 1 in (89), and therefore it becomes necessary to estimate an initial value of Q different from Q n at the first iteration. This value is initially estimated for m = 0, as the sum of the flow participating in the rest of vessels where SI n g < 1, as follows: allowing the transmission of a backward wave through the confluence.

Courant-Friedrichs-Lewy (CFL) Condition
The Godunov scheme in (15) is an explicit numerical scheme where wave celerities define the maximum allowable time step. Considering that the grid size ∆x is fixed, the time step ∆t is computed considering the speed of the faster wave, leading to the following limit: where ∆t min,i is the minimum time step among the time steps defined for each cell and ∆t min,JRP is the minimum time step among the time steps defined for each boundary cell involved in a JRP, where C k is the wave speed of the resulting non-linear wave, S k in (50) if the solution is a shock or λ 1g,k in (44) if the solution is a rarefaction. A CFL value equal to 1/2 ensures no interaction of waves from neighboring Riemann problems [26].

Results
In this section, many numerical tests are presented to show the capabilities of the JRP solver, comparing analytical and numerical solutions under extreme flow conditions on collapsible tubes, subsonic-supersonic transitions, sonic blockage conditions, and the transmission of backward waves from downstream vessels through the junction to upstream supersonic vessels with supersonic conditions. The analytical solutions provided by the JRP solver are compared with numerical results computed using the JRP at the junction in combination with the HLLS scheme. The solutions satisfy in all cases the entropy conditions for shocks or rarefactions.

JRP with Two Vessels
Prior to analyzing the numerical results for a JRP with multiple vessels or branches, in this section, we will first compare the exact solutions of a set of RPs with discontinuous geometrical and mechanical properties, with the analytical solutions provided by the JRP solver in (87) for the equivalent JRP with two vessels.
The exact solutions of the selected RPs involve non-linear waves and a stationary contact discontinuity generated by the variation of the mechanical and geometrical properties [25]. In some RPs, we also will plot the numerical solution of the RP using the Godunov updating scheme in (15), where the numerical fluxes in all edges will be computed exclusively using the HLLS scheme.
The equivalent JRP is defined assuming that two vessels with different geometrical and mechanical properties are connected through a junction. An analytical solution for the JRP is computed using the JRP solver in (87). We will also provide numerical solutions of the JRP using the Godunov updating scheme in (15). The JRP solver in (87) will provide the star values U k = [A k , Q k ] at the ending edges of the boundary cells. These start values will be used to compute fluxes according to equation (19). The HLLS scheme is used to compute the numerical fluxes at the inner edges of each vessel.
Initial conditions and mechanical and geometrical properties are presented in Table 1 for all the cases. In all cases reported in this work, density is ρ = 1 g/cm 3 Table 1. Two-vessel Junction configuration. Test case, number of vessel k, father or daughter vessel g k , reference diameter d o,k (mm), reference pulse wave velocity c o,k (cm s −1 ), reference pressure p o,k (mmHg), external pressure p e,k , initial speed index SI n = u n k /c n k , initial relative area α n k and transmural coefficients m and n. In test case 2V0, and in test cases 2V1 to 2V7, we check that the proposed JRP solver in (87) can recover solutions in the subsonic regime. The analytical solutions provided by the JRP solver, plotted in Figures 14-21 with a continuous line (-), are exactly equal to the exact solutions of the equivalent subsonic RPs in all these cases. In RpEnergyCase1 the initial solution is in exact equilibrium and the initial solution is preserved over time by the numerical solution, as expected. Test cases 2V1 to 2V7 are suitable combinations of subsonic compression/decompression backward/forward waves. The numerical solutions of the JRP also converge to the exact solutions as shown in Figures 14-21. Being the flow subsonic, in all these cases, conservation of total energy, p T , is observed at the junction.

Test Case
Test cases 2V8 and 2V9 involve supersonic flow. Figures 22-24 show the exact solutions and the numerical solutions when the problems are defined as a RP in a single vessel with discontinuous properties. Being the flow supersonic, the exact solution of RP RpCase8 ( Figure 22) is given by a contact wave at x = 0 followed by a BDW (backward compression wave) and an FDW (forward decompression wave). In RP 2V9, the exact solution, plotted in Figure 24, is a contact wave at x = 0 followed by an FCW (forward compression wave) and by an FDW (forward decompression wave). When solving both cases as a JRP using solver in (87), the analytical solution provided is different, since at most, one non-linear wave can be defined in each vessel of the confluence.
The analytical solution provided by system in (87) for JRP 2V8 ensures no change of the solution at the left vessel of the junction, as observed for the exact solution, but in the right vessel only a supersonic FDW is defined ( Figure 25). If we compare the analytical and the numerical solutions of the JRP, we observe supersonic flow conditions are kept in the left vessel, while in the right vessel, an FCW not defined in the analytical solution appears, followed by an FDW.
For test 2V9, the analytical solution of the JRP ensures no change of the solution at the left vessel of the junction, as in the exact solution, but in the right vessel only a supersonic FDW is defined. The numerical solution for the JRP is a supersonic flow entering the confluence, followed by an FCW (not defined analytical solution) followed by an FDW. It can be observed that, for both exact solutions of the RPs 2V8 and 2V9 total energy, p T , is preserved, while in JRPs 2V8 and 2V9, only mass conservation is preserved at the confluence in the analytical solution.           The pressure variation in JRP SubAtm1 leads to a subsonic left moving rarefaction wave and a progressive narrowing of the vessel area. A contact discontinuity appears at x = 0, leading to an abrupt expansion on the right side of the JRP, where sub-atmospheric pressure is imposed. In JRP SubAtm2 the suction on the left side of the JRP is increased and sonic conditions are forced across the contact discontinuity as shown in Figure 27. The moving rarefaction ends in the narrowest section now achieving sonic conditions. In both JRPs, SubAtm1 and SubAtm2, total energy is constant across the discontinuity. As sonic blockage has already been generated in JRP SubAtm2, the increase in the suction on the  The performance of the proposed JRP solver in (87) is checked in situations of vessel collapse in test cases Collapse1 to Collapse4 in [28]. Figures 29 and 30 show the exact solution of the RP Collapse1 in an artery. Supersonic initial flow conditions are imposed in all the vessel. In this RP, the solution remains constant over time on the left side of the RP. A stationary contact discontinuity is generated by the change in elasticity at x = 0, where a constant value of total energy is observed and then, on the right side of the RP, a supersonic BDW develops leading to an extreme reduction of the vessel area. The vessel area is recovered along a supersonic FDW wave. When the vessel collapses between the BDW and the FDW, the value of the PWV increases exponentially against small variations in area. Although the HLLS solver is capable of generating a suitable solution for flow or velocity in the collapsed region, the numerical solution for the SI differs from the exact one since the value reached for the area is not as low as the value reached by the value of the solution exact.
The solution changes when the JRP solver in (87) is applied to compute analytical solution for the test case Collapse1. The analytical solution does not change on the left side of the RP, but now, at the discontinuity, total energy does not remain constant. A supersonic FCW appears on the right side of the JRP, as only one wave can be defined at each branch of the confluence. The numerical solution to this JRP retains the exact solution in the left vessel, while in the right vessel, an FCW not defined in the analytical solution appears, followed by an FCW.  Test case Collapse2 has supersonic involves initial conditions on the left side and subsonic ones on the right. When solving this problem as a JRP, the analytical solution is equal to the exact solution of the equivalent RP in a single vessel, plotted in Figure 31. Numerical results can accurately reproduce the resulting supersonic BDW on the left and subsonic FDW on the right. JRP or equivalent RP Collapse3 is a vacuum problem where initial subsonic conditions lead to a supersonic BDW and a supersonic FDW, equal in both cases, as shown in Figure 32. When the solution of the JRM is computed numerically using the Godunov scheme, the results follow accurately the exact solution.
When solving RP Collapse4, initially at rest, the exact solution, shown in Figure 33, leads to a supersonic BDW on the left side of the RP. This rarefaction crosses the interface x = 0, ends on the right side and is followed by a supersonic FCW, as shown in Figure 33. When solving the problem as JRP the analytical solution leads to a subsonic BDW and a supersonic FCW on the left and right vessels respectively, as shown in Figure 34. Numerical solution of the JRP recovers accurately the analytical solution as shown in Figure 34.

JRP with Three Vessels
In this section, the analytical solution will be compared with the numerical solution for JRPs with three vessels. Initial conditions and mechanical and geometrical properties are presented in Table 2  JRPs 3V1 to 3V11 are RPs at junctions where three venous vessels are defined, one at the left side of the RP, and two at the right side. Initial conditions and mechanical and geometrical properties are equal at the right side of the RP. The analytical solution will be compared with the numerical solution for JRPs with three vessels.
In JRPs 3V1 to 3V3, vessel areas on the right side are larger than the vessel area on the left side. Initial conditions include zero velocity in all the domain. In absence of an external pressure, the variation in energy produced by the vessel deformation would generate a flow from the right vessels to the left vessel. However, in these cases the suction generated in the left side of the problem ensures a flow in the opposite direction, from the right vessel to the left vessels. Analytical solutions for the proposed JRP solver in (87) to test cases 3V1 to 3V3 are plotted in Figures 35-37 respectively, and they show the evolution of the JRP when suction is increased in the right vessels. While in test case 3V1 a subsonic rarefaction is generated in the left vessel, in test case 3V2 conditions of sonic blockage have achieved. Further suction in JRP 3V3 is unable to generate a large value of flow and solutions for JRPs 3V2 and 3V3 are equal. Numerical solutions of the JRPs using the Godunov updating scheme in (15) reproduce the solution accurately for these three cases.
In test case 3V4, the initial conditions for velocity in JRP 3V3 are modified and flow conditions close to the sonic regime are imposed in all vessels. The effect of the suction combined with the initial velocity generates again sonic blockage but now, in the right vessels, supersonic conditions are achieved along the FDW, as shown in Figure 38. Even the jump conditions in the FDW are subsonic, when the solution crosses by α PWV,min the minimum value of PWV is defined, and a maximum value of SI appears. Table 2. Three vessel confluence configuration. Test case, number of vessel k, father or daughter vessel g k , reference diameter d o,k (mm), reference pulse wave velocity c o,k (cm s −1 ), reference pressure p o,k (mmHg), external pressure p e,k , inital speed index SI n = u n k /c n k , initial relative area α n k and transmural coefficients m and n. In JRPs 3V5 and 3V6, test case 3V1 is modified and sonic conditions, SI = 1, and super-sonic conditions, SI = 1.2, are imposed, respectively. As shown in Figures 39 and 40, in both cases the initial condition does not change over time in the left vessel, resulting in a supersonic outflow boundary condition at the junction, while an FCW appears in the right vessels.

Test Case
In test case 3V7 the vessel area is equal to the reference vessel area in all vessels, α = 1, and subsonic flow conditions are imposed in all vessels (Figure 41). The resulting solution produces a subsonic-supersonic transition at the junction with the same value of total energy for the right vessels, with a subsonic BDW in the left vessel and a supersonic FDW in the right vessels. Numerical solutions do not converge to this solution, leading to a two-wave pattern in the right side vessels.
Subsonic and supersonic initial conditions in the left and right side of the JRP in test case 3V8 generate a subsonic BDW including sonic limitation and subsonic FCW. As plotted in Figure 42, the numerical solution generates a different level of energy in the vicinity of the junction if compared with the analytical solution. Therefore, it does not converge to the analytical solution and generates a BCW in the left vessel and an identical FCW for the right vessels.
A large pressure gradient is forced in test case 3V9 by imposing a large variation in the vessel deformation between the left vessel and the right vessels. The analytical solution in Figure 43 is a subsonic BDW in the left vessel, and a supersonic FCW in the right vessels, which share the same level of energy at the junction.
Initial supersonic conditions are imposed in all vessels in test case 3V10, as plotted in Figure 44. The analytical solution is a supersonic flow entering the confluence from the left vessel, and a supersonic FDW in the outlet vessels. Supersonic left vessel conditions remain unaltered over time and therefore left vessel only takes part in the mass conservation equation at the junction. Right vessels, experience the supersonic conditions at the junctions too, and share the same level of energy. The numerical solution reproduces the main features of the solution, although a small perturbation is visible at the vicinity of the junction. In test case 3V11, the flow moves in the opposite direction, from right to left, and again all vessels have initial supersonic conditions. Now initial conditions remain unaltered over time in the right vessels, and they only take part in mass conservation, whereas a supersonic BCW appears in the right vessel ( Figure 45).    x (cm) Figure 45. Section 3.2. Test 3V11. Comparison between analytical (-) and numerical solutions at t = 0.003 s. The solution from vessel 1 to 3 is: a supersonic BCW, a supersonic flow entering the confluence, and a supersonic flow entering the confluence, respectively. Total pressure is not preserved at x = 0.

JRP with Four Vessels
In the following test cases the analytical solution and the numerical solutions for RPs at junctions with four vessels are presented. In all test cases, the four vessels have the same mechanical and geometrical properties, listed in Tables 3-5 Table 3. Confluence configuration with four vessels. Test case, number of vessel k, father or daughter vessel g k , reference diameter d o,k (mm), reference pulse wave velocity c o,k (cm s −1 ), reference pressure p o,k (mmHg), external pressure p e,k , initial speed index SI n = u n k /c n k , initial relative area α n k and transmural coefficients m and n.  Table 4. Confluence configuration with four vessels. Test case, number of vessel k, father or daughter vessel g k , reference diameter d o,k (mm), reference pulse wave velocity c o,k (cm s −1 ), reference pressure p o,k (mmHg), external pressure p e,k , initial speed index SI n = u n k /c n k , initial relative area α n k and transmural coefficients m and n.  Table 5. Confluence configuration with four vessels. Test case, number of vessel k, father or daughter vessel g k , reference diameter d o,k (mm), reference pulse wave velocity c o,k (cm s −1 ), reference pressure p o,k (mmHg), external pressure p e,k , initial speed index SI n = u n k /c n k , initial relative area α n k and transmural coefficients m and n. In this JRP, 4V1, all vessels are initially at rest, each one with a different level of deformation, α (Figure 46). While in left vessels 1 and 2 the deformation parameter α > 1, in right vessels 3 and 4 we impose α < 1. The pressure discontinuity distribution is enlarged by the suction in right vessels 3 and 4 and generates a sonic BDW and a subsonic BDW for vessels 1 and 2, and subsonic FCWs for vessels 3 and 4. In vessel 1, sonic limitation is observed at the junction. Vessel 1 also exhibits the greatest level of energy at the junction, while vessels 2 to 4, provide the same level of energy at the confluence.

Test Case
In test case 4V2, the initial conditions are switched to generate the same exact wave pattern, but in the opposite direction, u ≤ 0, involving again sonic flow limitation in one vessel (Figure 47). The JRP solver presented in this work converges to the expected physically based solution.
In JRP 4V3 (Figure 48), the mechanical properties in 4V2 are changed by decreasing the PWV in all vessels. The same level of suction is prescribed for vessels 1 and 2. Now, the vessels are more easily deformable and, despite the level of suction has been reduced if compared with test case 4V2, sonic blockage appears in vessels 3 and 4 at the junction. Only vessels 1 and 2 share the same level of energy at the junction, where a subsonic BCW appears in each branch.
In test case 4V4 (Figure 49), the suction in vessel 2 is increased. Now vessels 1, 3 and 4 contribute to the flow mass towards vessel 2 ( Figure 49). As sonic blockage for vessels 1, 3 and 4 is produced, they exhibit different values of total energy at the junction. The solution for vessel 2 is a subsonic BCW.
In test case 4V5 (Figure 50), the large jump in the deformation parameter, α = 1.2 at the left and α = 0.2 at the right, generates a BDW for vessels 1 (ending in sonic limitation) and 2, and a supersonic FCW and a subsonic FCW for vessels 3 and 4, respectively. Subsonic vessels 2 and 4 retain the same value of energy at the junction shared by the supersonic FCW in inflow vessel 3.
When the jump in pressure in test case 4V5 is inverted, the solution, in test case 4V6 (Figure 51) is exactly equal, but in the opposite direction, u ≤ 0.
Test case 4V7 (Figure 52) is defined by reducing the minimum value of the deformation parameter to α = 0.1 in test case 4V6. The solution for vessels 3 and 4, on the right side of the junction, remains unaltered. The flow distribution changes on the right side generating a supersonic and a subsonic BCW for vessels 1 and 2, respectively. Except in vessel 3, where sonic limitation is produced, all vessels share the same value of energy being subsonic or supersonic.
In test case 4V8 (Figure 53), supersonic and subsonic conditions for inflow vessels 1 and 2 are prescribed. Rest conditions are imposed the right side of the JRP. The solution generates a subsonic FDW for vessel 4 and supersonic BDW for vessels 1 and 2, all of them connected by the same level of energy at the junction. Vessel 3, experiencing sonic blockage, participates in mass conservation but not in the energy state at the junction.
In test case 4V9 (Figure 54), the initial conditions are modified to generate the same solution as the one generated in test case 4V9, with u > 0.
In test case 4V10 (Figure 55), maximum deformation values are found at the right side of the JRP together with zero velocity. Minimum deformation values involving sonic and supersonic conditions are found at the left side of the JRP. Right vessels 3 and 4 develop a BDW, while vessels 1 and 2 experience the generation of subsonic BCW. The solution at x = 0 is subsonic for all branches, exhibiting the value of energy at the junction.
Maximum deformation values together with rest conditions are prescribed at the left side of the JRP in test case 4V11 (Figure 56). Supersonic and subsonic conditions with u < 0 are imposed on the right side. The jump in pressure generates a positive flow, restrained by the initial flow conditions in the right vessels. In right vessels 3 and 4 a subsonic FCW is developed, while in vessels 1 and 2 a subsonic BDW is generated. The solution at x = 0 is subsonic for all vessels.
Supersonic initial conditions are initially imposed on all vessels in test case 4V12 ( Figure 57) with u < 0, retaining maximum deformation values at the left and minimum deformation values at the right. The solution generates supersonic BDWs for vessels 1 and 2. Initial conditions in vessels 3 and 4 do not change over time. Even the flow is supersonic at x = 0 for vessels 1 and 2; they share the same level of p T at this point.
In test case 4V13 (Figure 58) ,supersonic initial conditions are initially imposed to define the same wave pattern as in test case 4V12, but in the opposite direction.
Supersonic initial conditions are imposed in outflow vessels 3 and 4 in test case 4V14 (Figure 59), where now the minimum vessel area has been defined, retaining their initial state over time, and therefore only contributing to the mass balance at the junction. Vessels 1 and 2 develop supersonic BDWs. They share the same level of energy at the junction despite ending at x = 0 in supersonic and subsonic conditions, respectively.
Subsonic initial conditions for vessel 3 and supersonic initial conditions for vessels 1, 2 and 4 are imposed in test case 4V15 (Figure 60). Vessels 1, 2, and 3 share the same value of p T at x = 0, and an FCW and two BDWs are generated for vessels 1, 2, and 3 respectively. Flow in vessel 4 does not change over time.
In 4V16 (Figure 61) supersonic initial conditions are initially imposed in vessels 1 and 2. While for vessel 2 the solution remains invariant over time, it can be seen how for vessel 1, a BCW is developed, allowing transmission of the information downstream the junction. Vessels 1, 3 and 4 shares the same level of energy at x = 0. In vessels 3 and 4, a subsonic FCW is observed.
Collapsed and partially collapsed initial vessel areas are imposed in test cases 4V17 to 4V21, 0.1 ≤ α ≤ 0.6. In test case 4V17 (Figure 62) supersonic conditions are imposed in vessels 2 to 4. While vessel 2 preserves its initial condition, the rest of the vessels accommodate their solutions to share a constant level of energy at the junction.
In test case 4V18 (Figure 63), subsonic/supersonic initial conditions are imposed in vessels 1,4/2,3, respectively. In the solution, while merging-flow vessel 2 preserves its initial condition, the rest of the vessels share a constant level of energy at the junction under conditions of subsonic flow.
Subsonic/supersonic initial conditions are imposed in vessels 1,4/2,3, respectively, in test case 4V19 (Figure 64), but with setting u < 0. Vessel 3 preserves its supersonic initial condition and does not participate in the energy level at the junction. The rest of the vessels share a constant level of energy at the junction under conditions of subsonic flow.
In test case 4V20 (Figure 65), the flow in vessels 1 and 4 is initially subsonic with u > 0, while the flow in vessels 2 and 3 is initially supersonic with u < 0. Flow in vessel 2 does not change over time. A supersonic BCW is generated in vessel 1 departing from subsonic conditions. A supersonic FDW and a subsonic FCW appear in vessels 3 and 4 respectively. Vessels 1, 3 and 4, even showing different conditions for the SI at x = 0, share the same level of energy at the junction.
The flow in vessels 1 and 4 is initially supersonic with u < 0, while the flow in vessels 2 and 3 is initially subsonic with u > 0, in test case 4V21 (Figure 66). Subsonic conditions for vessels 1 and 2 and supersonic conditions for vessel 3 are generated at the junction, under the same level of p T . The flow in Vessel 4 remains unaltered over time.
In test case 4V22 (Figure 67), supersonic initial conditions are imposed in all vessels, leading to a supersonic BCW and supersonic FDW for vessels 1 and 2 respectively, sharing the same level of p T at x = 0. Vessels 1 and 4 remain unaltered over time.
Zero velocity conditions are prescribed in vessel 2, while supersonic conditions are imposed in vessels 1,2 and 4, in test case 4V23 (Figure 68) . Solution in vessel 2, a BDW, generates sonic flow conditions at the junction, while supersonic initial flow conditions in vessel 1 remains unaltered over time. Both vessels 1 and 2 provide mass to vessels 3 and 4, where suction has been initially imposed. Flow in vessels 3 and 4 evolves from supersonic to subsonic conditions along an FDW and they share the same level of energy at the junction. It can also be clearly seen that for vessels 3 and 4, the SI does not show a monotonically decreasing variation between the initial and the solution region. The maximum value of SI appears when the vessel deformation becomes equal to α min during the expansion fan.
Test cases 4V24 to test case 4V27 are cases where collapsed areas are imposed in all vessels. Test case 4V23 is modified by imposing collapsed areas for vessels 3 and 4, leading to JRP 4V24 (Figure 69). Solutions for vessel 1 and 2, supersonic flow entering the confluence and a subsonic BDW including flow limitation, do not participate in the energy conservation equation. The initial conditions of the JRP generate an FDW for both vessels 3 and 4, and share the same level of p T at x = 0 under conditions of subsonic flow.
Initial conditions in test case 4V25 (Figure 70) ensure that all vessels are flow merging vessels. The solution is a supersonic flow entering the confluence for vessels 2 and 3 (unaltered over time), a subsonic BDW including sonic limitation for vessel 1 and a subsonic FCW for vessels 4, where a large jump in vessel area is observed. As a consequence, all vessels provide a different value of total pressure preserved at x = 0.
Again, all vessels are flow merging vessels in test case 4V26 (Figure 71). The solution is a supersonic flow entering the confluence for vessels 1 and 4 (unaltered over time), a subsonic BDW including sonic limitation for vessel 2 and a subsonic FCW for vessels 3, where now an expansion in vessel area is present. Total pressure is not preserved at x = 0.
Initial conditions in test case 4V27 ( Figure 72) generate a vacuum state in the neighborhood of the junction. All solutions are decompression waves ending in a subsonic state at x = 0, where all vessels share the same value of total pressure at x = 0.          x (cm) Figure 65. Section 3.3. Test 4V20. Comparison between analytical (-) and numerical solutions at t = 0.001 s. The solution of test case 4V20 from vessels 1 to 4 is: a supersonic BCW, a supersonic flow entering the confluence, a supersonic FDW, and a subsonic FCW, respectively. Vessels 1, 3, and 4 share the same value of p T at x = 0. x (cm) Figure 66. Section 3.3. Test 4V21. Comparison between analytical (-) and numerical solutions at t = 0.001 s. The solution of test case 4V21 from vessels 1 to 4 is: a supersonic BDW, a subsonic BCW, a supersonic FCW, and a supersonic flow entering the confluence, respectively. Vessels 1, 2, and 3 share the same value of p T at x = 0.    x (cm) Figure 72. Section 3.3. Test 4V27. Comparison between analytical (-) and numerical solutions at t = 0.001 s. The solution of test case 4V27 from vessels 1 to 4 is: a supersonic BDW, a subsonic BDW, a subsonic FDW, and a supersonic FDW, respectively. All vessels share the same value of p T at x = 0.

Conclusions
In the current work, a methodology for the solution of the Junction Riemann Problem for 1D hyperbolic balance laws in networks of elastic vessels has been proposed. The solver presented can deal with all possible transitions of flow, including flow limitation and supersonic flow conditions on elastic collapsible tubes. Solutions are based on rarefaction waves or shock waves that are completely consistent with the underlying hyperbolic system of conservation laws. The resulting JRP solver can be used in combination with any numerical scheme of an arbitrary order of accuracy.
The methodology has been illustrated in test cases with discontinuous variations of mechanical and geometrical properties of vessels and compared with exact solutions in junctions with two vessels, showing the expected limitations of the JRP solver presented in this work, where solutions in each branch of the confluence only contain at most one type of wave. The solutions satisfy, in all cases, the entropy conditions for shocks or rarefactions for both subsonic and supersonic conditions. The solver has also been tested in junctions with three and four vessels, considering merging flows, splitting flows or any other possible combination in a confluence of vessels under subsonic and supersonic conditions.
The solution method is based on an iterative search of the solution, and the computational details used in this work have been detailed. Convergence to suitable solutions has been assessed in a variety of cases. In all cases, analytical solutions of the proposed JRP solver and numerical solutions generated by combining the JRP methodology with the HLLS are presented. In all cases, convergence to the solution is guaranteed, including cases with collapsed states, subsonic-supersonic transitions, sonic blockage conditions, and the transmission of backward waves from downstream vessels through the junction to upstream supersonic vessels with supersonic conditions. Funding: The present work was partially funded by the Aragón Government through the Fondo Social Europeo.