EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6642 ISSN 1307-5543 – ejpam.com Published by New York Business Global Simulating Arterial Blood Flow Using Finite Volume Conservation Element and Solution Element Scheme Saqib Zia1,∗, Aqsa Qandeel1, Asad Rehman2, Sidrah Ahmed3, Munshoor Ahmed1, Wei Sin Koh4 1 Department of Mathematics, COMSATS University Islamabad, Park Road Chak Shehzad, Islamabad, Pakistan 2 Department of Information Technology, Pothohar Campus Gujar Khan, University of the Punjab, Pakistan 3 Department of Mathematics, Basic Sciences and Humanities, Sukkur IBA University, Sukkur, Pakistan 4 INTI International University, Persiaran Perdana BBN Putra Nilai, 71800 Nilai, Negeri Sembilan, Malaysia Abstract. This article investigates the numerical approximation of blood flow in arteries with spatially varying stiffness and resting cross-sectional areas. The underlying model is formulated to preserve steady-state solutions by accurately balancing non-zero fluxes with corresponding source terms. The study focuses on addressing the numerical challenges posed by the model’s complexity. To overcome these difficulties, a Conservation Element and Solution Element (CE/SE) scheme is employed, ensuring the preservation of the steady-state condition associated with zero discharge. The proposed scheme effectively minimizes numerical oscillations and dissipation, delivering im- proved accuracy and robustness compared to conventional methods. Its performance is evaluated through several benchmark cases from the existing literature, with results compared against those obtained using the Central Scheme. 2020 Mathematics Subject Classifications: 65N30, 74S05, 76M10 Key Words and Phrases: 1-D blood flow, CE/SE scheme, central scheme, KFVS scheme, non-conservative, hyperbolic systems, smart microgrids 1. Introduction Simulation of blood flow has numerous applications in medical science, driven by the pulsatile nature of blood and the necessity to accurately capture wave propagation within ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6642 Email addresses: saqibzia81@hotmail.com (S. Zia), aqsaqandeel18@gmail.com (A. Qandeel), assad013@gmail.com (A. Rehman), dr.sidrah@iba-suk.edu.pk (S. Ahmed), manshoor@comsats.edu.pk (M. Ahmed), weisin.koh@newinti.edu.my (W. S. Koh) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 2 of 18 arteries [1, 2]. The study of one-dimensional blood flow equations dates back to 1775, when Leonhard Euler first formulated them (as cited in [1]), highlighting the inherent challenges in their solution; see, for example, [2]. Since then, pulsatile wave flow has been extensively examined in seminal works, notably those of Pedley [3] and Lighthill [4]. Recent studies have made substantial progress in solving the governing systems of equations for arterial blood flow [5–17]. While advances in three-dimensional fluid–structure interaction models have been remarkable [18, 19], one-dimensional modeling remains indispensable particularly for analyzing complex arterial networks and venous systems [6–8, 20–23]. The shallow water equations (SWE), first formulated by Adhémar Barré de Saint- Venant in 1871 [24], have been extensively studied and applied across a wide range of fields. They are commonly used to model phenomena such as river flows [25], lake dynam- ics, dam-break scenarios, and tsunami propagation [26–28]. Incorporating source terms into the Saint-Venant equations, however, introduces numerical challenges, particularly in preserving steady-state solutions. To address this, Bermúdez and Vázquez [29, 30] pro- posed modifications to the Roe solver aimed at maintaining steady states. LeVeque [31] and Jin [32] developed exact and approximate Riemann solvers specifically designed for non-homogeneous problems, while other researchers [33, 34] proposed approaches based on solving the Riemann problem associated with the SWE. More recently, Delestre, Li, and collaborators [35–37] have applied finite volume schemes to solve the SWE for the simulation of blood flow in arteries. This study extends a CE/SE scheme [38] to address the SWE’s in a blood flow model for arteries. Unlike conventional methods, this scheme is specifically designed to overcome key limitations of existing approaches [38–41]. Its distinguishing characteristics include (i) a coherent treatment of time and space, (ii) the integration of conservation elements (CEs) and solution elements (SEs), and (iii) the ability to capture shocks without relying on Riemann solvers. The CE/SE scheme has been successfully applied in various domains, including unsteady flow dynamics [38–40], aeroacoustics [42], diffusion processes [41], and more [43–47]. The performance of the CE/SE scheme is assessed by comparing its results with those derived from the staggered central scheme [48]. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 3 of 18 Notation Symbol Description x Spatial coordinate along the vessel t Time Av(x, t) Cross-sectional area of the vessel A0(x) Resting cross-sectional area of the vessel Rv(x, t) Vessel radius R0 v Resting vessel radius Qv(x, t) Discharge (volume flow rate) ū(x, t) Mean blood velocity across the section u(x, t, r) Longitudinal velocity at radial location r ρ Blood density νb Blood viscosity Cf = 8πνb Friction coefficient k Arterial stiffness (elasticity parameter) W Conserved variable vector F (W ) Flux vector τ(W ) Source term vector J(W ) Jacobian matrix of flux λ1, λ2 Eigenvalues (wave speeds u± c) u Average flow velocity (Qv/Av) c Wave propagation speed Ω, V Conservation/solution element regions hi, fi Auxiliary flux/state vectors in CE/SE scheme ∆x,∆t Spatial and temporal step sizes α Limiter parameter in slope reconstruction N Number of computational grid cells L Length of the computational domain ϵ Perturbation amplitude for test problems ∆R Amplitude of radius variation (aneurysm cases) Φ Prescribed wave profile function c0 Characteristic wave speed Qamp Amplitude of imposed discharge wave ω Angular frequency of periodic inflow Tpulse Period of pulse wave kr, kj Real and imaginary components of wave number Table 1: List of symbols and notations used throughout the paper. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 4 of 18 2. Mathematical Equations The one-dimensional mathematical equations are derived from the Navier–Stokes (NS) equations by applying the long-wave approximation [49, 50]. In this formulation, trans- verse and longitudinal viscous effects are considered negligible. Incompressibility is ensured by assuming uniform pressure across each cross-section, while the flow variables vary only along the longitudinal coordinate (see Figure.1). The governing equations of the model are given by:  ∂Av ∂t + ∂Qv ∂x = 0 ∂Qv ∂t + ∂(Qv 2 Av ) ∂x + Av ρ ∂p ∂x = −Cf Qv Av . (1) Here, Av(x, t) = πRv 2 depicts the cross-sectional area of the vessel, Rv represents vessel radius,, Qv(x, t) = Av(x, t)ū(x, t) signifies the discharge, ū(x, t), is the mean flow velocity, ρ signifies blood density and Cf = 8πνb where νb defined as viscosity. ū(x, t) = ∫ Rv 0 2πu(t, x, r)rdr πRv 2 with u(t, x, r) signifies the longitudinal velocity. The 1D mathematical equations (1) can also be written as, By using the following relation, ū(x, t) = ∫ Rv 0 2πu(t, x, r)rdr πRv 2 with Av= πR0 2, A0 and k represent the resting cross-sectional area and arterial stiffness, respectively. The 1-D blood flow model (1) can be written in the following form, ∂Av ∂t + ∂Qv ∂x = 0 ∂Qv ∂t + ∂(Qv 2 Av + k 3ρ √ π Av 3/2) ∂x = kAv 2ρ √ π √ A0 ∂A0 ∂x − Cf Qv Av . (2) In compact form, ∂tW + ∂tF = τ (3) where W depicts conserved variable vectors, F = F (W ) is flux vector and τ = τ(W ) W = ( Av Qv ) , F(W) = ( Qv Qv 2 Av + k 3ρ √ π Av 3/2 ) (4) τ(W) = ( 0 kAv ρ √ π ∂x √ A0 − Cf Qv Av ) . (5) The system is hyperbolic for Av > 0. Therefore, we have ∂F (W ) ∂x = ( 0 1 k √ Av 2ρ √ π − Q2 v A2 v 2Qv Av ) · ∂ ∂x ( Av Qv ) = J(W ) · ∂W ∂x S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 5 of 18 The two eigenvalues of the Jacobian matrix J(W ) are real and distinct. λ1 = Qv Av −√ kRv 2ρ = u− c and λ2 = Qv Av + √ kRv 2ρ = u+ c. Figure 1: The left figure illustrates blood flow in which u, Rν , and R0 ν represent the velocity, radius, and radius at rest of the flow, respectively. The right figure depicts the shallow flow in which u, h, and z represent the velocity, height, and bottom topography of the flow, respectively. 3. Derivation of the scheme The CE/SE scheme serves as a powerful numerical tool for addressing various problems in continuum mechanics. Originally developed for the numerical computation of conser- vation laws governing fluid dynamics, this scheme stands out for its ability to minimize numerical dissipation and oscillations effectively. Here, 1-D CE/SE method [38] is derived for equations (1)–(5). Let (x0, x1) = (t, x) be the coordinates in 2-D Cartesian space (E2). The following form can be derived by applying the Gauss divergence theorem to Eq. (3).∮ Ω(ν) hi · dS = ∫ V τi dV (6) in which (a) hi := (Wi, Fi) T , i = 1, 2, i.e., for all i = 1, 2, Wi and Fi are the elements of hi in the specified t and x directions, (b) dS := dσ n, where dσ represents the area and n is the unit outward normal vector to the surface component on Ω(V ). We divide the whole domain into a union of conservation elements (CE’s) and solution elements (SE’s). Equation (6) is applied over the space-time region CE, which accommodates discontinuities in the flow variables. For each j, k = 0,±1 2 ,±1,±3 2 , . . ., there exists an SE corresponding to each (j, k), as depicted in Fig. 1 (dashed curve). For any point (x, t) ∈ SEk j , fi, hi, and Wi are estimated by f∗ i , h ∗ i , and W ∗ i , and these variables take the following form: W ∗ i (x, t) n j := (Wi) k j + (Wit) k j (t− tk) + (Wix) k j (x− xj), (7) and f∗ i (x, t) n j := (fi) k j + (fit) k j (t− tk) + (fix) k j (x− xj), (8) S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 6 of 18 Moreover, τ∗i (x, t) n j := τ ( W ∗ i ,W ∗ ix ) . (9) Figure 2: Staggered space-time grid Using the chain rule, we acquire (fix) k j = 4∑ n=1 (fi,n) k j (Wnx) k j , (10) (fit) k j = 4∑ n=1 (fi,n) k j (Wnt) k j , (11) where fi,n := ∂f i ∂Wn , i, n = 1, 2. (12) Note that (Wi) k j , (Wix) k j , and (Wit) k j are constants in SE(j, k). The numerical equiva- lents of these values are Wi, ∂Wi ∂x , and ∂Wi ∂t at (xj , t k). Since hi ≜ (Wi, fi), we define h∗i k j := ( W ∗ i k j , f ∗ i k j )T . (13) Moreover, for any (x, t) ∈ SEk j , W ∗ i , f ∗ i , and τ∗i satisfy Eq. (3) under the assumption of smoothness: ∂W ∗ i (x, t) ∂t + ∂f∗ i (x, t) ∂x = τ∗i (x, t). (14) S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 7 of 18 Using Eqs. (7) and (8), Eq. (14) becomes (Wit) k j = −(fix) k j + τi ( (Wi) k j , (Wix) k j ) . (15) We notice that (fi) k j depends on (Wi) k j , (fix) k j depends on (Wi) k j and (Wix) k j , (fit) k j depends on (Wi) k j and (Wit) k j , and Qi depends on (Wi) k j and (Wix) k j . See [46] for details. Here, fluxes are computed by the AUSMV method [51]. Using Eqs. (7)–(9) and (13), we get (Wi) k j = (Wi) k− 1 2 j± 1 2 ∓ ∆x 4 [ (Wix) k− 1 2 j± 1 2 + (Wix) k j ] ∓ ∆t ∆x [ (fi) k− 1 2 j± 1 2 − (fi) k j ] ∓ ∆t2 4∆x [ (fit) k− 1 2 j± 1 2 + (fit) k j ] ± [ (Ql) k− 1 2 j± 1 2 − (Ql) k j ] . (16) The summation gives (Wi) k j = 1 2 [ (Wi) k− 1 2 j− 1 2 + (Wi) k− 1 2 j+ 1 2 + (Si) k− 1 2 j− 1 2 − (Sil) k− 1 2 j+ 1 2 ] + (Ql) k− 1 2 j+ 1 2 − (Ql) k− 1 2 j− 1 2 , (17) where (Si) k− 1 2 j± 1 2 = ∆x 4 (Wix) k− 1 2 j± 1 2 + ∆t ∆x (Fi) k− 1 2 j± 1 2 + ∆t2 4∆x (fit) k− 1 2 j± 1 2 . (18) The following slope limiter is used to overcome oscillations: (Wix) k j = Ui ( (Wix−) k j , (Wix+) k j ; α ) , i = 1, . . . , 4, (19) where α ≥ 0 is a parameter, and Ui = x− |x+|α + x+ |x−|α |x+|α + |x−|α . (20) Furthermore, ( Wix+ )k j = (W ′ i ) k j+ 1 2 − (Wi) k j ∆x/2 , (21) ( Wix− )k j = (Wi) k j − (W ′ i ) k j− 1 2 ∆x/2 , (22) with ( W ′ i )k j± 1 2 = (Wi) k− 1 2 j± 1 2 + ∆t 2 (Wit) k− 1 2 j± 1 2 , i = 1, . . . , 4. (23) Combining Eqs. (17) and (19) completes the derivation of the 1-D CE/SE scheme. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 8 of 18 4. Numerical test problems This section provides a variety of numerical test cases from literature. Problem 1: The ideal tourniquet This case study is also addressed in [37, 52]. A tourniquet is applied and then imme- diately removed. The following initial conditions at t=0 are considered, Av(x) =  π(RL v ) 2 if x<0, π(RR v ) 2 otherwise. and Qv = 0. The computational domain [ − 4 100 , 4 100 ] m is divided into 200 grid points. Furthermore, we take the following values of the parameters: k = 1.0×107 kgm−2 s−2, ρ = 1.060×103 kgm−3, RL v = 5.0×10−3 m, RR v = 4.0×10−3 m. The results, shown in Figure 3, clearly indicate that the CE/SE scheme closely matches the reference solution. Furthermore, the L1-errors for all schemes are listed in Table 1, demonstrating that the CE/SE scheme achieves lower errors than the other methods. 4 4.2 4.4 4.6 4.8 5 ×10 -3 Radius at t=0.005 Reference CESE KFVS Central NT R a d iu s ( m ) -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 x(m) (a) Radius Rv -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 x(m) 0 1 2 3 4 5 6 7 8 D is ch ar ge ( m 2 /s ) ×10-5 Discharge at t=0.005 Reference KFVS WENO Central (b) Discharge Qv(x) Figure 3: The numerical solutions of the ideal tourniquet problem on a mesh with 200 cells at t = 0.005 seconds. (a) Radius Rv, (b) Discharge Qv(x). Problem 2: Wave Equation The purpose of this test study is to check the numerical scheme’s capability in sim- ulating perturbed steady solutions. The initial conditions are as follows [36, 37]: The cross-sectional area is defined as Av(x) =  π(R0 v) 2, if x ∈ [ 0, 4L10 ] ∪ [ 6L 10 , L ] , π(R0 v) 2 [ 1.0 + 1 200 sin ( π x− 4L 10 2L 10 )]2 , if x ∈ [ 4L 10 , 6L 10 ] , S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 9 of 18 Table 2: Problem 3: Comparison of L1-errors in the schemes. Radius N CE/SE KFVS Central(NT) 50 0.0141 0.2011 0.4510 100 0.0043 0.0181 0.3110 200 0.0022 0.0054 0.3139 400 0.0012 0.0028 0.2035 800 6.7120× 10−4 0.0023 0.1540 1600 2.4123× 10−4 6.4235× 10−4 0.0139 where the domain [0, 0.16] is divided into 200 points. Furthermore, we take: Q(x) = 0 m3 s−1, k = 1.0×108 kgm−2 s−2, ρ = 1.060×103 kgm−3, R0 v = 4.0×10−3 m, L = 1.6×10−1 m. The exact solutions are given by [37]: Rv(x, t) = R0 v + 1 400 [Φ(x− c0t) + Φ(x+ c0t)] , u(x, t) = − c0 400R0 v [−Φ(x− c0t) + Φ(x+ c0t)] . The numerical results are generated using the CE/SE scheme and the staggered cen- tral scheme and are presented in Figure 4. Results computed on 200 cells are shown at t = 0.002, 0.004, and 0.006 s. The outcomes from the proposed scheme closely match the reference solution. Problem 3: The man at eternal rest Here, we examine a non-flow scenario involving a variation in radiusR0(x), representing the condition of a deceased individual with an aneurysm. The initial data are: Rv = R0 v =  R̃, if x ∈ [0, x1] ∪ [x4, L], R̃+ 0.5∆R [ sin ( x−x1 x2−x1 π − π 2 ) + 1 ] , if x1 ≤ x ≤ x2, R̃+∆R, if x2 ≤ x ≤ x3, R̃+ 0.5∆R [ cos ( x−x3 x4−x3 π ) + 1 ] , if x3 ≤ x ≤ x4. The domain [0, L] m is divided into 200 points. At t = 0, u = 0, the values of other parameters are: k = 1.0×108 kgm−2 s−2, R̃ = 4.0×10−4 m, ρ = 1.060×103 kgm−3, ∆R = 1.0×10−3 m, L = 1.4×10−1 m, x1 = 1.0×10−2 m, x2 = 3.05×10−2 m, x3 = 4.95×10−2 m, x4 = 7.0×10−3 m. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 10 of 18 0 0.05 0.1 0.15 x(m) 4 4.005 4.01 4.015 4.02 ra di us ( m ) ×10-3 Profile of radius at time "0" sec Reference CESE Central 0 0.05 0.1 0.15 x(m) 4 4.002 4.004 4.006 4.008 4.01 ra di us ( m ) ×10-3 Profile of radius at time "0.002" sec Reference CESE Central 0 0.05 0.1 0.15 x(m) 4 4.002 4.004 4.006 4.008 4.01 ra di us ( m ) ×10-3 Profile of radius at time "0.004" sec Reference CESE Central 0 0.05 0.1 0.15 x(m) 4 4.002 4.004 4.006 4.008 4.01 ra di us ( m ) ×10-3 Profile of radius at time "0.006" sec Reference CESE Central Figure 4: Results of problem 2. From top left to bottom right: (a) Initial radius (Top left) (b) Radius at t = 0.002s (Top right) (c)Radius at t = 0.004s (Bottom left) (d) Radius at t = 0.006s (Bottom right) It is evident from Figure. 5 that our suggested scheme maintains the well-balanced prop- erty. Problem 4: Propagation of a pulse to an expansion Initially, we start by analyzing the case of a pulse in a section RR v that is expanding: AL v > AR v . The other related parameters are: ρ = 1060kgm−3, k = 1.0× 108 kgm−2 s−2, RL v = 5.0× 10−3 m, RR v = 4.0× 10−3 m, L = 1.6× 10−1 m, ∆R = 1.0× 10−3 m. On a small scale, we take the following decreasing shape: R0 v(x) =  RR v +∆R if x ∈ [0, x1], RR v + ∆R 2 [1.0 + cos ( x−x1 x2−x1 π ) ] if x ∈ [x1, x2], RR v else, with ϵ = 5.0× 10−3. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 11 of 18 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 x(m) -8 -6 -4 -2 0 2 4 6 8 ve lo ci ty ( m /s ) ×10-3 The man at eternal rest results at time "5" sec Reference CESE Figure 5: Representation of well-balanced CE/SE scheme. In Figure. 6, we compare the numerical outcomes to the reference solutions at t = 0.002s and t = 0.006s. The numerical solutions greatly approximate the referenced ones and are equivalent to those in [36]. Problem 5: Propagation of a pulse from an expansion In this case, all of parameters of problem 4 are considered except the initial radius, however, pulse propagation is considered from an expansion. Initially, the radius at t = 0 is given as, Rv(x) = { R0 v(x)[1.0 + ϵ sin ( 100 20Lπ(x− 15L 100 ) ) ] if x ∈ [15L100 , 35L 100 ], R0 v(x) else. Here, ϵ = 5.0×10−3. The numerical outcomes at t = 0.002s and t = 0.006s are depicted in Fig. 7. One observes that the numerical solutions closely predict the referenced ones and are analogous to those in [36]. Problem 6: Wave damping This test problem is similar to the Womersley [53] problem. An investigation of viscous damping term in the linearized momentum equation is carried out in this problem. A periodic signal with a constant segment at rest is considered at the inflow. The model is given by [37], { ∂tA+ ∂xQ = 0, ∂tQ+ ∂x( Q2 A + k 3ρ √ π A 3 2 ) = kA ρ √ π ∂x( √ A0)− Cf Q A . (24) In the above model, Cf = 8πv, the viscosity of the blood is represented by v. The domain S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 12 of 18 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 x(m) -5 0 5 10 15 20 R -R 0 ( m ) ×10 -6 Approximation of error at time "0" sec Reference CESE Central 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 x(m) -5 0 5 10 15 R -R 0 ( m ) ×10 -6 Approximation of error at time "0.002" sec Reference CESE Central 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 x(m) -5 0 5 10 15 R -R 0 ( m ) ×10 -6 Approximation of error at time "0.006" sec Reference CESE Central Figure 6: The numerical solutions of the propagation of a pulse to an expansion on a mesh with 100 cells: (a) Initial errors R − R0 at t = 0.0s (Top left) (b) Errors R − R0 at t = 0.002s (Top right) (c)Errors R − R0 at t = 0.004s (Bottom) [0, 3] is divided into 200 grid points. The radius Rv and discharge Qv at t = 0 are, Av(x) = π(R0 v) 2, Qv(x) = 0. The values of other related parameters are: ρ = 1060 kg.m−3, k = 1 × 108 Pa.m−1, R0 v = 4 × 10−3 m. For this example, the solution is obtained up to t = 25 s. A damping wave is obtained in the domain [36] Q(t, x) = { 0 if krx > ωt, Qamp sin(ωt− krx)e k iX if krx<ωt, (25) with kr = [ ω4 c40 + ( ωCf πR0 v 2c20 )2] 1 4 cos( 1 2 arctan(− Cf πR0 v 2ω )) , kj = [ ω4 c40 + ( ωCf πR0 v 2c20 )2] 1 4 sin( 1 2 arctan(− Cf πR0 v 2ω )) , ω = 2π/Tpulse = 2π/0.5s, S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 13 of 18 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 x(m) -0.5 0 0.5 1 1.5 2 2.5 3 R -R 0 ( m ) ×10 -5 Approximation of error at time "0.00" sec t=0 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 x(m) -5 0 5 10 15 20 R -R 0 ( m ) ×10 -6 Approximation of error at time "0.002" sec Reference CESE Central 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 x(m) -5 0 5 10 15 20 R -R 0 ( m ) ×10 -6 Approximation of error at time "0.006" sec Reference CESE Central Figure 7: The numerical solutions of the propagation of a pulse from an expansion on a mesh with 100 cells (a) Initial errors R−R0 at t = 0.0s (b) Errors R−R0 at t = 0.002s (c) Errors R−R0 at t = 0.004s. c0 = √√ A0k 2 √ πρ = √ R0 vk 2ρ . At x = 0 the incoming discharge is, Qb(0, t) = Qamp sin(ωt)m 3/s, where, Qamp = 3.45× 10−7m3.s−3 being the amplitude of inflowing discharge. The results are presented in Fig. 8. The CE/SE scheme agree well with the referenced solutions and those presented in [36]. 5. Conclusions In this study, a one-dimensional blood flow model for arterial circulation was numer- ically investigated using the conservation element and solution element(CE/SE) scheme. The inclusion of non-conservative terms in the governing equations complicates the balance between fluxes and source terms, introducing significant numerical challenges. The CE/SE scheme proved efficient in addressing these difficulties by accurately capturing sharp dis- continuities while maintaining local and global conservation. This capability makes it S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 14 of 18 0 0.5 1 1.5 2 2.5 3 x(m) -3 -2 -1 0 d is c h a rg e ( m 2 / s ) ×10 -7 Profile of discharge at time "25" sec with C f = 0 Reference CESE Central 0 0.5 1 1.5 2 2.5 3 x(m) -3 -2 -1 0 d is c h a rg e ( m 2 / s ) ×10 -7 Profile of discharge at time "25" sec with C f = 0.000022 Reference CESE Central 0 0.5 1 1.5 2 2.5 3 x(m) -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 d is c h a rg e ( m 2 / s ) ×10 -7 Profile of discharge at time "25" sec with C f = 0.000202 Reference CESE Central x[m] 0 0.5 1 1.5 2 2.5 3 D is ch ar ge [m 2 /s ] ×10-8 -14 -12 -10 -8 -6 -4 -2 0 2 Discharge at t = 25s with C f =0.005053 Reference CESE Central Figure 8: The numerical solutions of the wave damping on a mesh with 200 cells at t = 25 seconds. The damping of a discharge wave with Cf = 0 (upper left), Cf = 0.000022 (upper right), Cf = 0.000202 (lower left), andCf = 0.005053 (lower right) particularly suitable for wave phenomena in arterial systems that involve discontinuities or steep gradients, such as shock-like pulse waves and flow phase transitions. Numerical test cases demonstrated the robustness of the proposed approach under varying hemo- dynamic conditions, with the CE/SE scheme achieving superior accuracy and reduced oscillations compared to central and KFVS methods. Future work may extend this frame- work to more realistic physiological settings, including patient-specific arterial networks, viscoelastic wall models, and three-dimensional fluid–structure interactions. Such exten- sions will further enhance the applicability of the method in medical simulations and the study of cardiovascular pathologies. References [1] D. Xiu and S. J. Sherwin. Parametric uncertainty analysis of pulse wave propaga- tion in a model of a human arterial network. Journal of Computational Physics, 226(2):1385–1407, 2007. [2] K. H. Parker. A brief history of arterial wave mechanics. Medical & Biological Engineering & Computing, 47(2):111–118, 2009. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 15 of 18 [3] T. J. Pedley. The Fluid Mechanics of Large Blood Vessels. Cambridge University Press, Cambridge, 1980. [4] M. J. Lighthill and J. Lighthill. Waves in Fluids. Cambridge University Press, Cambridge, 2001. [5] A. Lucca, S. Busto, L. O. Müller, E. F. Toro, and M. Dumbser. A semi-implicit finite volume scheme for blood flow in elastic and viscoelastic vessels. Journal of Computational Physics, 495:112530, 2023. [6] H. G. Abdelwahed, M. A. Abdelrahman, A. F. Alsarhan, and K. Mohamed. Numerical simulating the blood flow model via nonhomogeneous riemann solver scheme. Partial Differential Equations in Applied Mathematics, 11:100845, 2024. [7] S. Jana and S. Kuila. Riemann solutions of two-layered blood flow model in arteries. International Journal of Non-Linear Mechanics, 156:104485, 2023. [8] Q. Zhang, W. Sheng, and T. Xiao. Riemann problem and godunov-type scheme for a two-layer blood flow model. Applied Mathematics Letters, 135:108437, 2023. [9] R. Mondal, M. Minhajul, and T. Raja Sekhar. Interactions between elementary waves and weak discontinuity in two-layer blood flow through artery. Physics of Fluids, 36(3), 2024. [10] S. Jana and S. Kuila. On the riemann problem and interaction of elementary waves for two-layered blood flow model through arteries. Mathematical Methods in the Applied Sciences, 47(1):27–46, 2024. [11] B. Ghitti, C. Berthon, M. H. Le, and E. F. Toro. A fully well-balanced scheme for the 1d blood flow equations with friction source term. Journal of Computational Physics, 421:109750, 2020. [12] G. I. Montecinos, A. Santacá, M. Celant, L. O. Müller, and E. F. Toro. Ader scheme with a simplified solver for the generalized riemann problem and an average eno reconstruction procedure. application to blood flow. Computers & Fluids, 248:105685, 2022. [13] S. Shagolshem, B. Bira, and K. V. Nagaraja. Analysis of shock wave propagation in two-layered blood flow model via lie symmetry. International Journal of Non-Linear Mechanics, page 104761, 2024. [14] W. Sheng and S. Xu. The riemann problem for the blood flow model with body force term. Advances in Mathematical Physics, 2024(1):2992241, 2024. [15] J. Alastruey, K. H. Parker, J. Peiró, and S. J. Sherwin. Lumped parameter out- flow models for 1-d blood flow simulations: Effect on pulse waves and parameter estimation. Communications in Computational Physics, 4(2):317–336, 2008. [16] M. Willemet, V. Lacroix, and E. Marchandise. Inlet boundary conditions for blood flow simulations in truncated arterial networks. Journal of Biomechanics, 44(5):897– 903, 2011. [17] E. Marchandise, N. Chevaugeon, and J. F. Remacle. Spatial and spectral super convergence of discontinuous galerkin method for hyperbolic problems. Journal of Computational and Applied Mathematics, 215(2):484–494, 2008. [18] M. A. Fernández, J. F. Gerbeau, and C. Grandmont. A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid. Interna- S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 16 of 18 tional Journal for Numerical Methods in Engineering, 69(4):794–821, 2007. [19] F. N. Van de Vosse, J. De Hart, C. H. G. A. Van Oijen, D. Bessems, T. W. M. Gunther, A. Segal, and F. P. T. Baaijens. Finite-element-based computational methods for cardiovascular fluid-structure interaction. Journal of Engineering Mathematics, 47(3- 4):335–368, 2003. [20] Y. Watanabe, T. Asada, P. Y. Lagre, M. Saito, Y. Ikenaga, and M. Matsukawa. One- dimensional model for propagation of a pressure wave in a model of the human arterial network: Comparison of theoretical and experimental. Journal of Biomechanical Engineering, 2011. [21] M. S. Olufsen, C. S. Peskin, W. Y. Kim, E. M. Pedersen, A. Nadim, and J. Larsen. Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions. Annals of Biomedical Engineering, 28(11):1281– 1299, 2000. [22] J. M. Fullana and S. Zaleski. A branched one-dimensional model of vessel networks. Journal of Fluid Mechanics, 621:183–204, 2009. [23] H. Tao, R. Husam, J. J. Dheyaa, K. S. Pradeep, F-H. A. Ali, H. A. H. Abbas, F. H. Dheyaa, M. Hiba, H. I. Ameer, and S. S. S. Narinderjit. Computational anal- ysis of pulsatile blood flow and heat transport dynamics in middle cerebral artery aneurysms in different body physiologies. International Journal of Modern Physics C, 36(6):2450236, 2025. [24] B. St Venant. Theory of the non-permanent movement of the waters with application to the floods of the rivers and to the introduction of the marees in their bed. Comptes Rendus de l’Académie des Sciences, 73:148–154, 1871. [25] J. Burguete, P. Garćıa-Navarro, and J. Murillo. Friction term discretization and limitation to preserve stability and conservation in the 1d shallow water model: Ap- plication to unsteady irrigation and river flow. International Journal for Numerical Methods in Fluids, 58(4):403–425, 2008. [26] A. Valiani, V. Caleffi, and A. Zanni. Finite volume scheme for 2d shallow-water equations. application to the malpasset dam-break. In The 4th CADAM Workshop, Zaragoza, pages 63–94, November 1999. [27] A. Valiani, V. Caleffi, and A. Zanni. Case study: Malpasset dam-break simulation using a two-dimensional finite volume method. Journal of Hydraulic Engineering, 128(5):460–472, 2002. [28] S. Popinet. Quadtree-adaptive tsunami modelling. Ocean Dynamics, 61(9):1261– 1285, 2011. [29] A. Bermudez and M. E. Vázquez. Upwind methods for hyperbolic conservation laws with source terms. Computers & Fluids, 23(8):1049–1071, 1994. [30] A. Bermúdez, A. Dervieux, J. A. Desideri, and M. E. Vázquez. Upwind schemes for the two-dimensional shallow water equations with variable depth using unstructured meshes. Computer Methods in Applied Mechanics and Engineering, 155(1-2):49–72, 1998. [31] R. J. LeVeque. Finite Volume Methods for Hyperbolic Problems, volume 31. Cam- bridge University Press, Cambridge, 2002. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 17 of 18 [32] S. Jin. A steady-state capturing method for hyperbolic systems with geometrical source terms. ESAIM: Mathematical Modelling and Numerical Analysis, 35(4):631– 645, 2001. [33] L. Gosse. A well-balanced flux-vector splitting scheme designed for hyperbolic systems of conservation laws with source terms. Computers & Mathematics with Applications, 39(9-10):135–159, 2000. [34] T. Gallout, J. M. Hérard, and N. Seguin. Some approximate godunov schemes to compute shallow-water equations with topography. Computers & Fluids, 32(4):479– 513, 2003. [35] O. Delestre, A. Ghigo, J. M. Fullana, and P. Y. Lagrée. A shallow water with variable pressure model for blood flow simulation. arXiv preprint arXiv:1509.01917, 2015. [36] O. Delestre and P. Y. Lagrée. A well-balanced finite volume scheme for blood flow simulation. International Journal for Numerical Methods in Fluids, 72(2):177–205, 2013. [37] G. Li, O. Delestre, and L. Yuan. Well-balanced discontinuous galerkin method and finite volume weno scheme based on hydrostatic reconstruction for blood flow model in arteries. International Journal for Numerical Methods in Fluids, 86(7):491–508, 2018. [38] S. C. Chang. The method of space-time conservation element and solution element: A new approach for solving the navier–stokes and euler equations. Journal of Com- putational Physics, 119(2):295–324, 1995. [39] S. C. Chang, X. Y. Wang, and C. Y. Chow. New developments in the method of space- time conservation element and solution element: Applications to two-dimensional time-marching problems, 1994. [40] S. C. Chang, X. Y. Wang, and C. Y. Chow. The space-time conservation element and solution element method: A new high-resolution and genuinely multidimensional paradigm for solving conservation laws. Journal of Computational Physics, 156(1):89– 136, 1999. [41] S. C. Chang, X. Y. Wang, and W. M. To. Application of the space-time conserva- tion element and solution element method to one-dimensional convection diffusion problems. Journal of Computational Physics, 165(1):189–215, 2000. [42] M. Liu, J. B. Wang, and K. Q. Wu. The direct aero-acoustics simulation of flow around a square cylinder using the ce/se scheme. Journal of Algorithms & Computational Technology, 1(4):525–538, 2007. [43] C. Y. Loh and K. B. M. Q. Zaman. Numerical investigation of transonic resonance with a convergent-divergent nozzle. AIAA Journal, 40(12):2393–2401, 2002. [44] C. Y. Loh, L. S. Hultgren, and S. C. Chang. Wave computation in compressible flow using space-time conservation element and solution element method. AIAA Journal, 39(5):794–801, 2001. [45] C. Loh, S. C. Chang, P. Jorgenson, and L. Hultgren. Noise computation of a shock- containing supersonic axisymmetric jet by the ce/se method. In 38th Aerospace Sci- ences Meeting and Exhibit, page 475, December 1999. [46] S. Qamar and S. Mudasser. On the application of a variant ce/se method for solving S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6642 18 of 18 two-dimensional ideal mhd equations. Applied Numerical Mathematics, 60(6):587– 606, 2010. [47] X. Y. Wang, C. L. Chen, and Y. Liu. The space-time ce/se method for solving maxwell’s equations in time-domain. In IEEE Antennas and Propagation Society International Symposium, volume 1, pages 164–167. IEEE, June 2002. [48] H. Nessyahu and E. Tadmor. Non-oscillatory central differencing for hyperbolic con- servation laws. Journal of Computational Physics, 87(2):408–463, 1990. [49] L. Formaggia, A. Quarteroni, and A. Veneziani. Cardiovascular mathematics: Mod- eling and simulation of the circulatory system. 2009. [50] M. S. Olufsen. A one-dimensional fluid dynamic model of the systemic arteries. Amer- ican Journal of Physiology-Heart and Circulatory Physiology, 276(1):H257–H268, 2000. [51] N. Wei, C. Xu, Y. Meng, G. Li, M. Xiao, and A. Liu. Numerical simulation of gas- liquid two-phase flow in wellbore based on drift flux model. Applied Mathematics and Computation, 338:175–191, 2018. [52] O. Delestre, C. Lucas, P. A. Ksinant, F. Darboux, C. Laguerre, T. N. T. Vo, and S. Cordier. Swashes: A compilation of shallow water analytic solutions for hydraulic and environmental studies. International Journal for Numerical Methods in Fluids, 72(3):269–300, 2013. [53] J. R. Womersley. XXIV. oscillatory motion of a viscous liquid in a thin-walled elastic tube: The linear approximation for long waves. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 46(373):199–221, 1955.