An HDG-type discretization of 1D embedded blood flow in complex networks with Riemann-coupled junctions

We present a hybridizable discontinuous Galerkin (HDG) method for one-dimensional blood flow in arterial networks with Windkessel boundary conditions. The scheme uses a Harten–Lax–van Leer (HLL) numerical flux and implicit time integration. After spatial discretization, the complete network—including cell states, hybrid traces, junction conditions, and dynamic $\mathrm{RCR}$ terminal pressures—is written as a residual index-one differential–algebraic equation (DAE) and integrated with variable-step, variable-order backward differentiation formulas. The implicit formulation removes the CFL stability restriction of explicit schemes and allows variable-order, variable-step solvers. The resulting single global residual and Jacobian for the network and its terminal dynamics also provide a natural interface for monolithic coupling with external nonlinear models. We establish local solvability and a semi-discrete energy balance, and show numerically optimal convergence rates for both solution variables in the $L^2$ norm and the broken $H^1$ seminorm. The solver is validated against established benchmarks, including an aortic bifurcation and a $37$-artery network, and is further applied to a cerebral arterial network. Mass conservation is observed in a synthetically generated anatomical liver network of $8\,\mathrm{k}$ vessels, and strong-scaling studies are performed on synthetic networks of $50\,\mathrm{k}$ and $100\,\mathrm{k}$ vessels.