Wave propagation in stenotic vessels; theoretical analysis and   comparison between 3D and 1D fluid-structure-interaction models

George Papadakis; Jean Raspaud

arXiv:1905.04173·physics.flu-dyn·July 24, 2019

Wave propagation in stenotic vessels; theoretical analysis and comparison between 3D and 1D fluid-structure-interaction models

George Papadakis, Jean Raspaud

PDF

TL;DR

This paper develops an analytical 1D model for wave propagation in stenotic vessels, accurately approximating 3D FSI results and enabling fast, large-scale simulations for cardiovascular studies.

Contribution

The authors introduce a simple, analytically solvable 1D model using only two geometric parameters that closely matches complex 3D FSI simulations for stenotic vessels.

Findings

01

Excellent agreement with 3D FSI volume flow rate (<1% error)

02

Accurate pressure distribution prediction with small error

03

Model is adaptable for arbitrary vessel shapes and large-scale simulations

Abstract

Using analytical expressions for the pressure and velocity waveforms in tapered vessels, we construct a linear 1D model for wave propagation in stenotic vessels in the frequency domain. We demonstrate that using only two parameters to approximate the exact geometry of the constriction (length and degree of stenosis), we can construct a model that can be solved analytically and can approximate with excellent accuracy the response of the original vessel for a wide range of physiologically relevant frequencies. We then proceed to compare the 1D results with full 3D FSI results from the literature for parameters corresponding to an idealized stenotic carotid artery. We find excellent matching with the volume flow rate over the cardiac cycle (less than $1%$ error). Using results from DNS simulations to parametrize the velocity profile in the stenotic region, we manage to predict also the…

Tables2

Table 1. Table 1 : Fluid and wall properties.

Fluid Density $ρ_{f} (k g / m^{3})$	$1000$
Wall Density $ρ_{s} (k g / m^{3})$	$1000$
Modulus of elasticity $E (P a)$	$10^{6}$
Poisson ratio $ν (-)$	0.3
Wall Thickness $h (m)$	$0.5 \times 10^{- 3}$

Table 2. Table 2 : Fluid and membrane parameters of the 3D FSI model (data taken from Figueroa et al. ( 2006 ) ).

Modulus of Elasticity ( $d y n / c m^{2}$ )	$4.07 \times 10^{6}$
Fluid density ( $g / c m^{3}$ )	$1.06$
Membrane density ( $g / c m^{3}$ )	$1.0$
Membrane thickness ( $c m$ )	$0.03$
Poisson ratio	$0.5$
Inlet/Outlet vessel diameter ( $c m$ )	$0.6$
Throat diameter ( $c m$ )	$0.3$
Length of stenotic region ( $c m$ )	$1.0$

Equations63

\frac{\partial}{\partial s} (s^{2} u_{s}) + (f s^{3}) \frac{\partial p}{\partial t} = 0

\frac{\partial}{\partial s} (s^{2} u_{s}) + (f s^{3}) \frac{\partial p}{\partial t} = 0

\frac{\partial u _{s}}{\partial t} = - \frac{1}{ρ} \frac{\partial p}{\partial s}

f = \frac{( 1 - ν ^{2} ) tan ^{2} α sin α}{E h ( 1 - cos α )}

f = \frac{( 1 - ν ^{2} ) tan ^{2} α sin α}{E h ( 1 - cos α )}

\frac{\partial}{\partial t} [p Q_{s}] + [0 2 π (1 - cos α) \frac{s ^{2}}{ρ} \frac{1}{2 π ( 1 - c o s α ) f s ^{3}} 0] \frac{\partial}{\partial s} [p Q_{s}] = [00]

\frac{\partial}{\partial t} [p Q_{s}] + [0 2 π (1 - cos α) \frac{s ^{2}}{ρ} \frac{1}{2 π ( 1 - c o s α ) f s ^{3}} 0] \frac{\partial}{\partial s} [p Q_{s}] = [00]

R (s) = [- Z (s) 1 Z (s) 1]

R (s) = [- Z (s) 1 Z (s) 1]

Y (s) = \frac{1}{Z ( s )} = 2 π (1 - cos α) \frac{f}{ρ} s^{+ 5/2}

Y (s) = \frac{1}{Z ( s )} = 2 π (1 - cos α) \frac{f}{ρ} s^{+ 5/2}

w_{1} (s) = \frac{1}{2 λ ρ} (p + λ ρ u_{s}), w_{2} (s) = \frac{1}{2 λ ρ} (- p + λ ρ u_{s})

w_{1} (s) = \frac{1}{2 λ ρ} (p + λ ρ u_{s}), w_{2} (s) = \frac{1}{2 λ ρ} (- p + λ ρ u_{s})

p (x, t) = P (x, ω) e^{iω t}, u_{s} (x, t) = U_{s} (x, ω) e^{iω t},

p (x, t) = P (x, ω) e^{iω t}, u_{s} (x, t) = U_{s} (x, ω) e^{iω t},

P (x, ω)

P (x, ω)

U_{s} (x, ω)

Q_{s} (x, ω) = - iY (s) \frac{1}{s} (A J_{\frac{4}{3}} (z) + B Y_{\frac{4}{3}} (z))

Q_{s} (x, ω) = - iY (s) \frac{1}{s} (A J_{\frac{4}{3}} (z) + B Y_{\frac{4}{3}} (z))

P (s_{2}^{i}) = P (s_{1}^{i + 1}) \Rightarrow \frac{1}{s _{1}^{i}} (A_{i} J_{\frac{1}{3}} (z_{2}^{i}) + B_{i} Y_{\frac{1}{3}} (z_{2}^{i})) = \frac{1}{s _{1}^{i + 1}} (A_{i + 1} J_{\frac{1}{3}} (z_{1}^{i + 1}) + B_{i + 1} Y_{\frac{1}{3}} (z_{1}^{i + 1}))

P (s_{2}^{i}) = P (s_{1}^{i + 1}) \Rightarrow \frac{1}{s _{1}^{i}} (A_{i} J_{\frac{1}{3}} (z_{2}^{i}) + B_{i} Y_{\frac{1}{3}} (z_{2}^{i})) = \frac{1}{s _{1}^{i + 1}} (A_{i + 1} J_{\frac{1}{3}} (z_{1}^{i + 1}) + B_{i + 1} Y_{\frac{1}{3}} (z_{1}^{i + 1}))

U_{s} (s_{2}^{i}) = U_{s} (s_{1}^{i + 1}) \Rightarrow (\frac{f ^{i}}{ρ})^{\frac{1}{2}} (A_{i} J_{\frac{4}{3}} (z_{2}^{i}) + B_{i} Y_{\frac{4}{3}} (z_{2}^{i})) = (\frac{f ^{i + 1}}{ρ})^{\frac{1}{2}} (A_{i + 1} J_{\frac{4}{3}} (z_{1}^{i + 1}) + B_{i + 1} Y_{\frac{4}{3}} (z_{1}^{i + 1}))

U_{s} (s_{2}^{i}) = U_{s} (s_{1}^{i + 1}) \Rightarrow (\frac{f ^{i}}{ρ})^{\frac{1}{2}} (A_{i} J_{\frac{4}{3}} (z_{2}^{i}) + B_{i} Y_{\frac{4}{3}} (z_{2}^{i})) = (\frac{f ^{i + 1}}{ρ})^{\frac{1}{2}} (A_{i + 1} J_{\frac{4}{3}} (z_{1}^{i + 1}) + B_{i + 1} Y_{\frac{4}{3}} (z_{1}^{i + 1}))

U_{s} (s_{2}^{i_{t h r}}) = - U_{s} (s_{1}^{i_{t h r} + 1}) \Rightarrow (\frac{f ^{i_{t h r}}}{ρ})^{\frac{1}{2}} (A_{i_{t h r}} J_{\frac{4}{3}} (z_{2}^{i_{t h r}}) + B_{i_{t h r}} Y_{\frac{4}{3}} (z_{2}^{i_{t h r}})) = - (\frac{f ^{i_{t h r} + 1}}{ρ})^{\frac{1}{2}} (A_{i_{t h r} + 1} J_{\frac{4}{3}} (z_{1}^{i_{t h r} + 1}) + B_{i_{t h r} + 1} Y_{\frac{4}{3}} (z_{1}^{i_{t h r} + 1}))

U_{s} (s_{2}^{i_{t h r}}) = - U_{s} (s_{1}^{i_{t h r} + 1}) \Rightarrow (\frac{f ^{i_{t h r}}}{ρ})^{\frac{1}{2}} (A_{i_{t h r}} J_{\frac{4}{3}} (z_{2}^{i_{t h r}}) + B_{i_{t h r}} Y_{\frac{4}{3}} (z_{2}^{i_{t h r}})) = - (\frac{f ^{i_{t h r} + 1}}{ρ})^{\frac{1}{2}} (A_{i_{t h r} + 1} J_{\frac{4}{3}} (z_{1}^{i_{t h r} + 1}) + B_{i_{t h r} + 1} Y_{\frac{4}{3}} (z_{1}^{i_{t h r} + 1}))

\frac{1}{s _{1}^{1}} (A_{1} J_{\frac{1}{3}} (z_{1}^{1}) + B_{1} Y_{\frac{1}{3}} (z_{1}^{1})) = P_{0}

\frac{1}{s _{1}^{1}} (A_{1} J_{\frac{1}{3}} (z_{1}^{1}) + B_{1} Y_{\frac{1}{3}} (z_{1}^{1})) = P_{0}

w_{2} (s_{1}^{N}) = 0 \Rightarrow - P (s_{1}^{N}) + λ (s_{1}^{N}) ρ U_{s} (s_{1}^{N}) = 0 \Rightarrow - \frac{1}{s _{1}^{N}} (A_{N} J_{\frac{1}{3}} (z_{1}^{N}) + B_{N} Y_{\frac{1}{3}} (z_{1}^{N})) + (ρ f^{N} s_{1}^{N})^{- \frac{1}{2}} ρ [\frac{1}{i} (\frac{f ^{N}}{ρ})^{\frac{1}{2}} (A_{N} J_{\frac{4}{3}} (z_{1}^{N}) + B_{N} Y_{\frac{4}{3}} (z_{1}^{N}))] = 0

w_{2} (s_{1}^{N}) = 0 \Rightarrow - P (s_{1}^{N}) + λ (s_{1}^{N}) ρ U_{s} (s_{1}^{N}) = 0 \Rightarrow - \frac{1}{s _{1}^{N}} (A_{N} J_{\frac{1}{3}} (z_{1}^{N}) + B_{N} Y_{\frac{1}{3}} (z_{1}^{N})) + (ρ f^{N} s_{1}^{N})^{- \frac{1}{2}} ρ [\frac{1}{i} (\frac{f ^{N}}{ρ})^{\frac{1}{2}} (A_{N} J_{\frac{4}{3}} (z_{1}^{N}) + B_{N} Y_{\frac{4}{3}} (z_{1}^{N}))] = 0

p_{4} (x, t) = p_{4}^{f} e^{iω (t - \frac{x - x _{3}}{c _{4}})} + p_{4}^{b} e^{iω (t + \frac{x - x _{3}}{c _{4}})} = p_{4}^{f} e^{iω (t - \frac{x - x _{3}}{c _{4}})}

p_{4} (x, t) = p_{4}^{f} e^{iω (t - \frac{x - x _{3}}{c _{4}})} + p_{4}^{b} e^{iω (t + \frac{x - x _{3}}{c _{4}})} = p_{4}^{f} e^{iω (t - \frac{x - x _{3}}{c _{4}})}

\displaystyle Q_{4}\left(x,t\right)=Y_{4}\left(p_{4}^{f}e^{i\omega(t-\frac{x-x_{3}}{c_{4}})}-p_{4}^{b}e^{i\omega(t+\frac{x-x_{3}}{c_{4}})}\right)=Y_{4}p_{4}^{f}e^{i\omega(t-\frac{x-x_{3}}{c_{4}})}\

p_{3} (s, t) = s^{- \frac{1}{2}} (A_{3} J_{\frac{1}{3}} (z) + B_{3} Y_{\frac{1}{3}} (z)) e^{iω t}

p_{3} (s, t) = s^{- \frac{1}{2}} (A_{3} J_{\frac{1}{3}} (z) + B_{3} Y_{\frac{1}{3}} (z)) e^{iω t}

Q_{s}_{3} (s, t) = - iY (s) s^{- \frac{1}{2}} (A_{3} J_{\frac{4}{3}} (z) + B_{3} Y_{\frac{4}{3}} (z)) e^{iω t}

p_{4}^{f} = s_{3 L}^{- \frac{1}{2}} (A_{3} J_{\frac{1}{3}} (z_{3 L}) + B_{3} Y_{\frac{1}{3}} (z_{3 L}))

p_{4}^{f} = s_{3 L}^{- \frac{1}{2}} (A_{3} J_{\frac{1}{3}} (z_{3 L}) + B_{3} Y_{\frac{1}{3}} (z_{3 L}))

Y_{3} p_{4}^{f} = - iY (s_{3 L}) s_{3}^{- \frac{1}{2}} (A_{3} J_{\frac{4}{3}} (z_{3 L}) + B_{3} Y_{\frac{4}{3}} (z_{3 L}))

\frac{B _{3}}{A _{3}} = - \frac{i J _{\frac{4}{3}} ( z _{3 L} ) + J _{\frac{1}{3}} ( z _{3 L} )}{Y _{\frac{1}{3}} ( z _{3 L} ) + i Y _{\frac{4}{3}} ( z _{3 L} )}

\frac{B _{3}}{A _{3}} = - \frac{i J _{\frac{4}{3}} ( z _{3 L} ) + J _{\frac{1}{3}} ( z _{3 L} )}{Y _{\frac{1}{3}} ( z _{3 L} ) + i Y _{\frac{4}{3}} ( z _{3 L} )}

Y_{e f f} (x_{2}) = - iY (s_{2 R}) \frac{J _{\frac{4}{3}} ( z _{2 R} ) + \frac{B _{3}}{A _{3}} Y _{\frac{4}{3}} ( z _{2 R} )}{J _{\frac{1}{3}} ( z _{2 R} ) + \frac{B _{3}}{A _{3}} Y _{\frac{1}{3}} ( z _{2 R} )}

Y_{e f f} (x_{2}) = - iY (s_{2 R}) \frac{J _{\frac{4}{3}} ( z _{2 R} ) + \frac{B _{3}}{A _{3}} Y _{\frac{4}{3}} ( z _{2 R} )}{J _{\frac{1}{3}} ( z _{2 R} ) + \frac{B _{3}}{A _{3}} Y _{\frac{1}{3}} ( z _{2 R} )}

p_{2} (s, t) = s^{- \frac{1}{2}} (A_{2} J_{\frac{1}{3}} (z) + B_{2} Y_{\frac{1}{3}} (z)) e^{iω t}

p_{2} (s, t) = s^{- \frac{1}{2}} (A_{2} J_{\frac{1}{3}} (z) + B_{2} Y_{\frac{1}{3}} (z)) e^{iω t}

Q_{s}_{2} (s, t) = - iY (s) s^{- \frac{1}{2}} (A_{2} J_{\frac{4}{3}} (z) + B_{2} Y_{\frac{4}{3}} (z)) e^{iω t}

s_{2 L}^{- \frac{1}{2}} (A_{2} J_{\frac{1}{3}} (z_{2 L}) + B_{2} Y_{\frac{1}{3}} (z_{2 L})) e^{iω t} = p_{3} (s_{2 R}, t)

s_{2 L}^{- \frac{1}{2}} (A_{2} J_{\frac{1}{3}} (z_{2 L}) + B_{2} Y_{\frac{1}{3}} (z_{2 L})) e^{iω t} = p_{3} (s_{2 R}, t)

iY (s_{2 L}) s_{2 L}^{- \frac{1}{2}} (A_{2} J_{\frac{4}{3}} (z_{2 L}) + B_{2} Y_{\frac{4}{3}} (z_{2 L})) = Q_{s}_{3} (s_{2 R}, t) = Y_{e f f} (x_{2}) p_{3} (s_{2 R}, t)

\frac{B _{2}}{A _{2}} = \frac{i J _{\frac{4}{3}} ( z _{2 L} ) - \frac{Y _{e f f} ( x _{2} )}{Y ( s _{2 L} )} J _{\frac{1}{3}} ( z _{2 L} )}{\frac{Y _{e f f} ( x _{2 R} )}{Y ( s _{2 L} )} Y _{\frac{1}{3}} ( z _{2 L} ) - i Y _{\frac{4}{3}} ( z _{2 L} )}

\frac{B _{2}}{A _{2}} = \frac{i J _{\frac{4}{3}} ( z _{2 L} ) - \frac{Y _{e f f} ( x _{2} )}{Y ( s _{2 L} )} J _{\frac{1}{3}} ( z _{2 L} )}{\frac{Y _{e f f} ( x _{2 R} )}{Y ( s _{2 L} )} Y _{\frac{1}{3}} ( z _{2 L} ) - i Y _{\frac{4}{3}} ( z _{2 L} )}

Y_{e f f} (x_{1 R}) = iY (s_{1 R}) \frac{J _{\frac{4}{3}} ( z _{1 R} ) + \frac{B _{2}}{A _{2}} Y _{\frac{4}{3}} ( z _{1 R} )}{J _{\frac{1}{3}} ( z _{1 R} ) + \frac{B _{2}}{A _{2}} Y _{\frac{1}{3}} ( z _{1 R} )}

Y_{e f f} (x_{1 R}) = iY (s_{1 R}) \frac{J _{\frac{4}{3}} ( z _{1 R} ) + \frac{B _{2}}{A _{2}} Y _{\frac{4}{3}} ( z _{1 R} )}{J _{\frac{1}{3}} ( z _{1 R} ) + \frac{B _{2}}{A _{2}} Y _{\frac{1}{3}} ( z _{1 R} )}

\frac{p _{1}^{b}}{p _{1}^{f}} = \frac{1 - \frac{Y _{e f f} ( x _{1 R} )}{Y _{1}}}{1 + \frac{Y _{e f f} ( x _{1 R} )}{Y _{1}}} e^{- i 2 ω \frac{x _{1} - x _{0}}{c _{1}}}

\frac{p _{1}^{b}}{p _{1}^{f}} = \frac{1 - \frac{Y _{e f f} ( x _{1 R} )}{Y _{1}}}{1 + \frac{Y _{e f f} ( x _{1 R} )}{Y _{1}}} e^{- i 2 ω \frac{x _{1} - x _{0}}{c _{1}}}

Y_{e f f} (x_{0}) = Y_{1} \frac{1 - \frac{p _{1}^{b}}{p _{1}^{f}}}{1 + \frac{p _{1}^{b}}{p _{1}^{f}}}

Y_{e f f} (x_{0}) = Y_{1} \frac{1 - \frac{p _{1}^{b}}{p _{1}^{f}}}{1 + \frac{p _{1}^{b}}{p _{1}^{f}}}

p_{1} (x, t) = P (x_{0}) [cos (ω \frac{x - x _{0}}{c _{1}}) - i \frac{Y _{e f f} ( x _{0} )}{Y _{1}} sin (ω \frac{x - x _{0}}{c _{1}})] e^{iω t}

p_{1} (x, t) = P (x_{0}) [cos (ω \frac{x - x _{0}}{c _{1}}) - i \frac{Y _{e f f} ( x _{0} )}{Y _{1}} sin (ω \frac{x - x _{0}}{c _{1}})] e^{iω t}

\frac{\partial U _{x}}{\partial t} + U_{x} \frac{\partial U _{x}}{\partial x} = - \frac{1}{ρ} \frac{\partial P}{\partial x} + \frac{f}{ρ A}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Wave propagation in stenotic vessels; theoretical analysis and comparison between 3D and 1D fluid-structure-interaction models

G. Papadakis

Department of Aeronautics, Imperial College London, SW7 2AZ, UK

[email protected]

J. Raspaud

Abstract

Using analytical expressions for the pressure and velocity waveforms in tapered vessels, we construct a linear 1D model for wave propagation in stenotic vessels in the frequency domain. We demonstrate that using only two parameters to approximate the exact geometry of the constriction (length and degree of stenosis), we can construct a model that can be solved analytically and can approximate with excellent accuracy the response of the original vessel for a wide range of physiologically relevant frequencies. We then proceed to compare the 1D results with full 3D FSI results from the literature for parameters corresponding to an idealized stenotic carotid artery. We find excellent matching with the volume flow rare over the cardiac cycle (less than $1\%$ error). Using results from DNS simulations to parametrize the velocity profile in the stenotic region, we manage to predict also the pressure distribution with small error (a few percentage points). The method proposed in the paper can be used to approximate vessels of arbitrary shape profile and can be extended to cover the whole cardiovascular tree. Recursive expressions make the solution very fast and open the possibility of carrying out sensitivity and uncertainty quantification studies that require thousands (or even millions) of simulations with minimal cost.

keywords:

Stenotic vessels , Pulse wave propagation , 1-D blood flow modelling , comparison between 3D and 1D models.

††journal: Journal of Fluids and Structures

1 Introduction

Stenosis is the narrowing of arteries possibly due to inflammation or plaque buildup (known as atherosclerosis). In the coronary artery, stenosis can cause angina and is an established cause of heart failure; in the renal artery can result in kidney failure. The flow through a constricted vessel is considered as difficult to simulate, mainly because of complexity arising from unsteady flow separation and transition.

Early research into stenotic flows dates back to 1970’s, where experimental in vitro studies for steady and unsteady inflow conditions were conducted in idealized models by Young and Tsai (1973a, b). For a Reynolds number between 100 and 5000 and steady flow conditions the authors encountered three flow regimes: for low Reynolds numbers, the flow was laminar without separation, for larger numbers a zone of back flow was observed, and further increasing the Reynolds number, led to highly turbulent flow and no localized discrete region of reverse flow was discernible. They also found the flow characteristics are strongly affected by the percentage area-reduction of the stenosis. Numerous experimental studies have followed (Cassanova and Giddens (1978), Ahmed and Giddens (1983), Ojha et al. (1989), Karri and Vlachos (2010), Pielhop et al. (2012), Choi et al. (2017) to cite but a few).

A large number of 2D or 3D numerical simulations have been also conducted. Only 3D simulations can capture the loss of axisymmetry and transition (Sherwin and Blackburn (2005), Varghese et al. (2007a, b), Pal et al. (2014)). The aforementioned numerical studies ignore the wall distensibility, i.e. they assume a rigid wall. Three dimensional models that account for fluid-structure-interaction are those of Gerbeau et al. (2005) and Figueroa et al. (2006). Fully resolved 3D simulations are not suitable for modelling the whole arterial network. For this purpose, simpler 1D models have been successfully employed (Mynard and Nithiarasu (2008), Alastruey et al. (2011), Sherwin et al. (2003)). These models are formulated in physical, i.e. space-time, variables. More details on modeling studies (as well as experiments) on steady and unsteady flows in arteries are provided in the review articles of Berger and Jou (2000), Ku (1997), Taylor and Draney (2004) and van de Vosse and Stergiopulos (2011).

The central aim of the paper is to derive, apply and validate a linear 1D model for a stenotic vessel in terms of space-frequency (as opposed to space-time) variables. The model is based on analytical solutions of pressure and velocity waveforms in elastic, tapered vessels that were derived by Papadakis (2011). These solutions were found to match very well with 2D fluid-structure interaction (FSI) results. We present here the application of these solutions to model wave propagation in stenotic vessels and we validate the results against 3D FSI solutions.

The most significant advantage of the proposed analytical, frequency domain, solution method is that it is very fast. When the equations are formulated in space-time variables, the large speed of blood pulse waves and short lengths of vessels, restricts the time step to very small values, to the order of $10^{-4}s$ (Du et al. (2016)) or even less. Furthermore, one needs to evolve the solution for a sufficient number of cardiac cycles (each with period around $0.8s-1s$ ) to obtain an asymptotic periodic result. The combination of small time step and long integration times leads to an appreciable overall computing time.

For some applications, many simulations of the 1D model (thousands or even millions) may be necessary. Uncertainty quantification (see Chen et al. (2013)) and sensitivity analysis due to lack of accurate information on vessel properties, both require the application of the 1D model multiple times. To give an example, Huberts et al. (2013) investigated the effect of uncertainty of 73 parameters of the pulse propagation model. Their objective was to determine the most influential of these parameters when the model is applied to support decision making in arteriovenous fistula (AVF) surgery for hemodialysis. Their computational framework involved a variance-based method and Monte Carlo simulations. In total, they had to perform more than 2 million simulations, each taking an average 20s.

This is exactly where frequency domain methods offer a distinct advantage. For example, Flores et al. (2016) presented and validated a linear 1D formulation of blood flow in networks of flexible vessels for which a full analytical solution exists in the frequency domain. They call this the generalized Darcy’s elastic model (GDEM). The model is based on a linear relationship of the volume flow rate and pressure gradient in the frequency domain, known as Darcy’s law. Using appropriate compatibility conditions between successive vessels and at bifurcation points, they derive a linear system of equations with unknowns the values of pressure at the nodes of the cardiovascular tree (a node is defined as the mathematical point where two or three vessels are connected). Our approach is along the lines of Flores et al. (2016), and can be extended to cover the whole cardiovascular tree.

The paper is organized as follows. In section 2, we explain how the analytical solutions derived for tapered vessels can be used to model wave propagation in stenotic vessels, and in section 3 we apply the method to compute the frequency response of a vessel. In section 4, we simplify the form of the stenotic region and show that we can derive analytically solutions that match very closely with the results of the full model for physiologically relevant frequencies (section 5). In section 6, we compare the results of the developed 1D model against benchmark 3D solutions, and we conclude in section 7.

2 Mathematical formulation of 1D wave propagation model in a stenotic vessel

In this section we propose an analytical 1D model for wave propagation in a stenotic vessel. The profile of the vessel is approximated by a series of $N$ tapered segments (or elements), each of length $\Delta x$ , as shown in figure 1(a) below. The $i$ -th element is defined by the boundary faces $i-1$ and $i$ , $i=1...N$ , refer to figure 1(b).

The governing equations for each element take the form (for details refer to Papadakis (2011)):

[TABLE]

The equations are formulated in a local spherical reference frame, and the radius $s$ is measured from the origin of the coordinate system $O_{(i)}$ , as shown in figure 1(b). For the converging element $i$ shown in the figure, $s_{2}^{i}\leq s\leq s_{1}^{i}$ , where the subscripts $1$ and $2$ refer to flow inlet and outlet respectively. Velocity $u_{s}$ (positive in the $s$ direction) and pressure $p$ are cross-section averaged and $\rho$ is the density of the fluid. The area of the cross section is $A(s)=2\pi(1-\cos\alpha)s^{2}$ where $\alpha$ is the (local) cone angle. The flow is along the positive $x$ direction, so in the converging part of the vessel, $s$ and $x$ are in opposite directions, while in the expanding part they are at the same direction. Parameter $f$ in equation (1a) is a constant, and is defined as:

[TABLE]

It groups together the geometrical and wall parameters of the local segment (the tapering angle $\alpha$ , the modulus of elasticity $E$ , the thickness of the vessel wall $h$ , and the Poisson ratio $\nu$ ).

Introducing the volume flux $Q_{s}=A(s)u_{s}=2\pi(1-\cos\alpha)s^{2}u_{s}$ , system (1) can be written in matrix form as

[TABLE]

The eigenvalues of the coefficient matrix are the propagation velocities $\lambda_{1,2}=\pm\left(fs\rho\right)^{-1/2}$ . It can be easily shown that the matrix of eigenvectors $R(s)$ can be written as

[TABLE]

where $Z(s)=\dfrac{1}{2\pi(1-\cos\alpha)}\sqrt{\dfrac{\rho}{f}}s^{-5/2}$ is called resistance. The reciprocal

[TABLE]

is called admittance. The characteristic variables $w_{1},w_{2}$ propagating in the positive and negative $s$ directions respectively are,

[TABLE]

These variables are not invariants along a characteristic as in straight vessels; they are instead fully coupled because of the continuous reflection from the tapering wall.

It was shown in Papadakis (2011) that if the system of equations (1) is formulated in the frequency domain, i.e. we assume a solution of the form,

[TABLE]

then closed form analytic expressions can be obtained for the complex amplitudes $P\left(x,\omega\right)$ and $U_{s}\left(x,\omega\right)$ , as follows

[TABLE]

and after some algebra

[TABLE]

In the above expressions, ${\ J}_{\frac{1}{3}}\left(z\right),{\ Y}_{\frac{4}{3}}\left(z\right)$ are Bessel functions of order $\frac{1}{3}$ and $\frac{4}{3}$ respectively, with real argument $z=\frac{2}{3}\omega{(\rho f)}^{\frac{1}{2}}s^{\frac{3}{2}}$ . Constants $A$ and $B$ are determined by the boundary conditions imposed at the two ends of each segment, locations $s=s_{1}^{i}$ and $s=s_{2}^{i}$ (refer to figure 1(b)). These expressions incorporate analytically the change in amplitude of the forward and backward reflecting waves due to tapering. In the work of Flores et al. (2016), a linear relation between the volume flow rate and pressure gradient is applied (Darcy’s law). Such a relation is valid for cylindrical vessels only; for tapered vessels such an expression cannot be derived.

In total there are 2N unknown constants, the set of coefficients $A$ and $B$ for each segment (denoted below as $A_{i}$ and $B_{i}$ respectively for the $i-th$ segment). To evaluate these coefficients we apply

$2(N-1)$ compatibility conditions for velocity and pressure at the interfaces between two consecutive segments, and

2.

$2$ boundary conditions at the inlet and outlet of the whole vessel.

Figure 2(a) below shows two consecutive elements, $i$ and $i+1$ , that share the $i$ -th face and taper in the same direction (they both converge for increasing $x$ ). If the common face is located at the throat of the vessel, then the elements taper in opposite directions (i.e. one will converge and the other will expand), as shown in figure 2(b).

As mentioned earlier, radii $s_{1}^{i}$ and $s_{2}^{i}$ denote the inlet (upstream) and outlet (downstream) edges of the segment $i$ respectively. The compatibility condition for pressure between segments $i$ and $i+1$ leads to

[TABLE]

Similarly for velocity we get the condition

[TABLE]

For the interface located exactly at the throat, i.e. for $i=i_{thr}$ (refer to figure 2(b)), a change of sign is necessary for velocity i.e.

[TABLE]

At the inlet we assume that pressure is prescribed i.e.

[TABLE]

where $P_{0}=1$ , i.e. unit pressure excitation. For the outlet boundary, we assume that there is no incoming wave (non-reflecting condition) i.e. the characteristic variable $w_{2}(s_{1}^{N})=0$ in (6), which leads to

[TABLE]

Other boundary conditions can be easily implemented at the inlet and outlet. For example, in section 6 we apply a prescribed velocity at the inlet and impedance boundary condition at the outlet. The resulting system of equations for $A_{i},B_{i}$ has a lower block diagonal structure, with block size $2\times 2$ . Although in practice it was solved numerically, the structure of the system allows a recursive computation, with successive elimination of unknowns starting from the $N-th$ segment and proceeding backwards to the first. This is an important observation and will be evoked later in section 4 where a closed form analytical expression is derived for a simplified stenotic vessel profile. We note that in the work of Flores et al. (2016), the linear system does not have this convenient matrix structure, and can be solved only numerically.

3 Frequency response of a stenotic vessel

We apply the formulation of the previous section to compute the frequency response of a stenotic vessel that consists of three parts: a straight duct of length and radius $R=1cm$ , a symmetric converging-diverging nozzle with a total length $4R$ and a straight duct with length $5R$ . A sketch is shown in figure 3 and table 1 summarises the fluid and wall properties. These properties are valid for human vessels such as the splenic veins, Ünsal et al. (2006).

The vessel is discretised using $N=100$ tapered segments, each $1mm$ long. The converging-diverging part is approximated with $40$ elements (these values lead to grid independent results). The two straight parts are approximated as tapered vessels with a very shallow angle ( $\alpha=0.001^{\circ}$ ) to avoid a singularity in the evaluation of parameter $f$ (refer to (2)). The linear system was formulated and solved with $LU$ decomposition in MATLAB.

The degree of stenosis, $S$ , is defined as the normalized cross-sectional area reduction between the throat and the inlet i.e. $S=1-\left(\frac{R_{min}}{R}\right)^{2}$ . Three cases are considered with $S$ equal to $36\%$ , $51\%$ and $64\%$ . For each case, the vessel profile was obtained by fitting a $6th$ order polynomial through the following 7 points with (axial, radial) coordinates equal to $(0,R)$ , $(R,R)$ , $(2R,R_{mid})$ , $(3R,R_{min})$ , $(4R,R_{mid})$ , $(5R,R)$ and $(6R,R)$ , where $R_{mid}=0.5(R+R_{min})$ . A straight part with shallow angle was then added to extend the vessel to a total length of $10R$ .

The magnitude of the pressure and velocity response is plotted against the angular frequency $\omega$ at 4 spatial locations in figures 4 and 5 respectively. The superscript “ $FR$ ” in $P^{FR}$ and $U^{FR}$ in the vertical axes denotes ”Frequency Response” for unit pressure excitation at the inlet, i.e. $P_{0}=1Pa$ .

In order to interpret the results, we consider the angular frequency corresponding to a wavelength equal to a characteristic length scale of the vessel. For the case considered, we take this length scale to be $5R$ ; this is a dynamically important scale because there is no reflection in the straight part between $5R-10R$ . The wave propagation velocities at the straight parts and the throat are equal to $c_{str}=\sqrt{\frac{Eh}{(1-\nu^{2})\rho(2R)}}=5.24m/s$ and $c_{thr}=\sqrt{\frac{Eh}{(1-\nu^{2})\rho(2R_{min})}}=5.86m/s,6.26m/s$ and $6.76m/s$ for $S=36\%,51\%$ and $64\%$ respectively. Considering the average velocities $c_{ave}=5.55m/s,5.75m/s$ and $6.0m/s$ for the 3 values of $S$ and length of $5R$ we get the reference angular frequencies $\omega_{ref}=\frac{2\pi c_{ave}}{5R}=697rad/s,722rad/s$ and $753rad/s$ .

Both pressure and velocity plots show resonance peaks in a narrow frequency range, $[436,442]rad/s$ . In this range, the wavelength of the inlet excitation is of the same order as the length scale of the vessel and resonance occurs. This is due to the superposition of the forward and backward moving waves from the converging and diverging parts. Backward waves also reflect at the inlet and propagate forward. At every location therefore there is a multitude of waves that are superimposed. If they are in phase, a constructive interference occurs resulting in a significant amplification of the pressure, i.e. resonance. As expected, the amplification increases with the stenosis degree $S$ . It also depends on the location along the vessel; the converging part of the vessel amplifies the pressure wave more compared to the diverging part. Indeed, the maximum amplification appears at the middle of the converging part and it is equal to 7.5 for $S=64\%$ . Regarding velocity, the maximum amplification is located at the throat, as it can be expected (refer to figure 5).

The frequency range at which resonance occurs is however much larger that the physiological range, say between $[0,50]rad/s$ . For such low frequencies, the wavelength is much larger that the characteristic length scale of the vessel and the interaction is much weaker.

In order to obtain a deeper physical insight on the interaction between the forward and reflected waves in a stenotic vessel, in the following section we derive an analytical solution by simplifying the profile shape. We then compare the results of the simplified model to the full model in section 5.

4 Analytical solution for a simplified stenotic vessel.

We consider the simplified profile shape shown in figure 6 below. We have retained only the two most important geometric features, the length of the stenotic section and the degree of stenosis, $S$ (marked with solid blue line in the figure). We replace the detailed profile with two tapered and two straight segments, and therefore the whole vessel consists of 4 compartments. As will be shown below, this simplification allows us to derive a recursive analytic solution by combining expressions valid for straight and tapered sections. For the case examined in the previous section, the characteristic interface locations $x_{0},x_{1},x_{2},x_{3},x_{4}$ take the values $0,R,3R,5R$ and $10R$ respectively.

For the analytic solution, we need also the spherical radii of the inlet and outlet of each tapered compartment. These are denoted by “ $s$ ” followed by the index of the interface location and the letter “L” or “R”, depending on whether they were computed from the left or the right tapered section, see figure 6. For example, $s_{2L}$ is the spherical radius evaluated at interface 2 (between compartments 2 and 3) and is referring to the outlet of the left compartment (i.e. 2). The corresponding argument of the Bessel function is denoted as $z_{2L}=\frac{2}{3}\omega{(\rho f_{2})}^{\frac{1}{2}}s_{2L}^{\frac{3}{2}}$ .

Starting with compartment 4, we can write (see Lighthill (1975)):

[TABLE]

where $p_{4}^{f}$ and $p_{4}^{b}$ are the (constant) amplitudes of the forward and backward propagating waves, $c_{4}$ is the wave propagation velocity and $Y_{4}=\frac{A_{4}}{\rho c_{4}}$ is the admittance of compartment 4. Due to the non-relective boundary condition at $x=x_{4}$ , amplitude $p_{4}^{b}=0$ .

For the tapering compartment 3, we apply the analytic solutions (8a) and (9),

[TABLE]

Applying continuity of flow rate and pressure at the interface $x_{3}$ between compartments 3 and 4 we get

[TABLE]

Assuming that the admittance at the interface $x_{3}$ is continuous, i.e. $Y_{3}=Y(s_{3L})$ , and eliminating $p_{4}^{f}$ we get

[TABLE]

Strictly speaking, we expect a weak discontinuity of $Y$ at $x_{3}$ because of the slope change. We chose to ignore it here in order to simplify the resulting analytical expressions.

The effective admittance at $x_{2}$ is defined as $Y_{eff}\left(x_{2}\right)=\dfrac{Q_{s3}(s_{2R},t)}{p_{3}(s_{2R},t)}$ and is equal to

[TABLE]

with the ratio $\frac{B_{3}}{A_{3}}$ given by (18). We proceed in the same way for compartment 2, where

[TABLE]

and apply continuity at the interface $x_{2}$

[TABLE]

Notice that in the left hand side of (21b) the sign was changed from $(-)$ to $(+)$ . This is because $x_{2}$ is located at the throat of the stenotic vessel and the spherical radii in compartments 2 and 3 grow in opposite directions (this equation is the exact analogue of (12)). We take as positive the +x direction. After some algebra we get

[TABLE]

where the ratio $\frac{Y_{eff}\left(x_{2R}\right)}{Y\left(s_{2L}\right)}$ can be computed from equation (19) above, because again we assume continuity, i.e. $Y\left(s_{2L}\right)=Y\left(s_{2R}\right)$ . Now we can compute the effective admittance at $x_{1}$ , $Y_{eff}\left(x_{1}\right)=\dfrac{-{Q_{s}}_{2}(s_{1R},t)}{p_{2}(s_{1R},t)}$ (notice again the sign change in ${Q_{s}}_{2}$ ) as

[TABLE]

We apply the same procedure for compartment 1. The ratio $\frac{p_{1}^{b}}{p_{1}^{f}}$ is found to be

[TABLE]

and the effective admittance at $x_{0}$ is

[TABLE]

We are now able to complete the analytical solution. At the boundary $x_{0}$ , $p_{1}^{f}+p_{1}^{b}=P(x_{0})=1$ (assuming unit pressure excitation), and taking into account (25), we find that the pressure at an arbitrary point $x$ in the 1st compartment is given by,

[TABLE]

Similarly we can compute the pressure and volume flow rate at every point along the vessel.

Summarising, in order to evaluate the analytic solution for every $\omega$ , we start form the last compartment and proceed backwards. More specifically, we first compute $\frac{B_{3}}{A_{3}}$ from (18), and then $\frac{Y_{eff}(x_{2R})}{Y(s_{2R})}$ from (19), $\frac{B_{2}}{A_{2}}$ from (22), $\frac{Y_{eff}(x_{1R})}{Y(s_{1R})}$ from (23), $\frac{p_{1}^{b}}{p_{1}^{f}}$ from (24) and $\frac{Y_{e}ff(x_{0})}{Y_{1}}$ from (25) in that order. To the best of our knowledge, these analytic solutions have not appeared in the literature before.

The fact that a recursive solution can be obtained is directly related to the lower block diagonal structure of the system matrix, as mentioned in section 2. It is important to note that the same recursive approach can be easily extended to handle the whole cardiovascular tree. Appropriate compatibility conditions (pressure continuity and mass conservation) are applied at the tree junctions, and these allow the computation of the local admittances. More details can be found in Lighthill (1975).

5 Comparison of the analytic solution with the numerical solution of the full model.

The frequency response of the full model is compared with the analytic solution obtained in the previous section in figures 7 and 8 for pressure and velocity respectively.

It is very interesting to see that the simplified model can capture very accurately the frequency response of the full model for both pressure and velocity at all locations in a frequency range less than about $~{}400rad/s$ . This range is well within the physiological regime. The resonance frequency is also very well predicted, but the value is slightly underestimated.

These results indicate that only two parameters, namely the length and degree of stenosis, suffice to capture the behaviour of the vessel at the relevant frequencies. The exact details of the geometry affect the magnitude of the resonance peak and the response at large frequencies only.

We can easily explain this behaviour. When the frequency of the input signal is small, the wavelength is large and so the response is determined by the large-scale features of the stenotic vessel (degree of stenosis and length). On the other hand, larger frequencies that correspond to shorter wavelengths, are affected by the finer details of vessel, for example the exact stenotic profile. Of course, this is within the context of linear, inviscid theory. The effect of viscosity will be examined in the next section.

6 Comparison between 1D and 3D models

To assess the performance of the 1D model, we compare against the results of a 3D FSI model of an idealised stenotic carotid artery (see Figueroa et al. (2006)). The 3D model employs a membrane formulation for the vessel wall, which is considered as a linear elastic material. The membrane formulation, without transverse shear, was also employed to derive equations (1) of the 1D model. The fluid and structure equations are solved using the coupled momentum (or transpiration) method on a fixed computational domain. Discretisation is based on the finite element methods for both media.

The physical parameters are provided in table 2. The throat is located at $2.5cm$ from the inlet and the total length of the vessel is $12.6cm$ . The degree of stenosis is 75% and the ratio of stenotic length to inlet diameter, $L/D=1.67$ .

At the inlet, the volume flow rate is imposed (shown as the blue line in figure 10(a)). During the cardiac cycle, the Reynolds number based on the inlet diameter, $D$ , and the cross-section average velocity, $u_{m}(t)$ , varies between $213$ and $767$ , the time-average being $365$ (the peak-to-mean ratio is $3.6$ ). The reduced velocity $U_{red}=\overline{u_{m}(t)}T/D=42$ , where the overbar denotes time- averaging and $T=1.1s$ is the period of the cardiac cycle. Under these conditions, the flow downstream of the throat was found to be transitional and non-axisymmetric. This is consistent with the results of stability analysis of a stenotic vessel with the same $S$ but slightly larger $L/D$ ( $2$ instead of $1.67$ ) of Sherwin and Blackburn (2005). For steady inlet flow, the authors observed transitional flow for Reynolds larger than $688$ (for our case, the $Re$ exceeds this value for part of the cardiac cycle). The transition was found to be due a sub-critical bifurcation that leads to loss of axisymmetry and a weak Coanda-type deflection of the jet emanating from the throat towards the walls (similar behavior is shown in Figure 7 of Figueroa et al. (2006)). Sherwin and Blackburn (2005) also considered transition due to pulsating inlet flow, but their reduced velocities are much smaller (up to 7.5) compared to the present case. It is therefore more appropriate to compare with the steady results, as the reduced velocity is quite large.

In order to apply the 1D model, the flow rate is decomposed in Fourier series, with fundamental period $T$ . Most of the energy is included in the fundamental mode and the first 10 harmonics. The response for each harmonic was computed in the frequency domain using the approach of section 2 and the result was converted back to the time domain using inverse Fourier transform. The steady component was treated separately; the time-average volume flow rate was imposed at the inlet and the velocity at each location was computed using the continuity equation.

An impedance boundary condition was imposed at the outlet. More specifically, the volume and pressure signals at the outlet were written in Fourier series and the complex impedance $Z(\omega_{i})=P(\omega_{i})/Q(\omega_{i})$ was computed for each frequency $\omega_{i}$ . The imposed boundary condition was then $P(\omega_{i})=Z(\omega_{i})AU_{s}(\omega_{i})$ , where $A=\pi D^{2}/4$ is the outlet area. This approach reproduces the exact boundary condition of the 3D FSI model and eliminates any uncertainties in the comparison that can arise from inconsistent outlet conditions.

The fluctuating volume flow rate and pressure at the inlet and outlet of the domain during the entire cardiac cycle are shown in figure 9. There is a phase difference between the inlet and outlet due to the finite propagation velocity. This is not constant during systole because of the variation of the propagation velocity due to the presence of stenosis. Note also that during diastole, the volume flow rate at the outlet is larger compared to the inlet. This is because volume has accumulated inside the flexible vessel during systole, and released during diastole. This behaviour can only be captured using a model that accounts for vessel distensibility.

Adding the steady component of velocity to the fluctuating part, we obtain the total volume flow rate, which is compared to the 3D results in figure 10(a). At the outlet the values predicted by the 1D model match very closely with the 3D values during the entire periodic cycle. Indeed, the error between the 1D and 3D models (dashed and solid red lines respectively) is less than $1\%$ . Note that phase difference between inlet and outlet mentioned above is captured very accurately.

The computation of pressure is more involved. In the momentum equation of the 1D model we have ignored the non-linear and the viscous terms. Both are expected to be significant in the present case. In order to take them into account approximately, we write the 1D momentum equation as,

[TABLE]

where $f(x,t)$ is the frictional force per unit length and we have assumed that the profile shape factor is equal to $1$ (for details refer to Alastruey et al. (2012)). Assuming the velocity profile,

[TABLE]

that satisfies the non-slip condition $u(x,R,t)=0$ at the wall, we find that $f(x,t)=-2(\zeta+2)\mu\pi U_{x}(x,t)$ , where $\mu$ is the dynamic viscosity. Parameter $\zeta$ determines the shape of the profile. The standard parabolic profile corresponds to $\zeta=2$ . Increasing $\zeta$ makes the profile flatter around the centerline and increases shear close to the wall. Using the values of $U_{x}(x,t)$ from the 1D inviscid model (that were shown to match well the 3D results) and substituting in (27), we obtain the pressure gradient $\frac{\partial P}{\partial x}$ . In order to compute $P(x,t)$ , a boundary condition is needed and a value of $\zeta$ . Using the time-average value of $P$ from the 3D results at the outlet and adding the fluctuating pressure shown in figure 9(b) we obtain the required boundary condition. At the outlet the volume flow rate is well predicted, and since a boundary condition that matches the impedance between the 3D and 1D model was applied, we expect that the pressure profile will be also well predicted.

This is indeed the case, as can be clearly seen in figure 10(b) (red solid and dashed lines overlap). This confirms that our implementation is correct. The results in this figure were obtained with $\zeta=2$ for the whole domain. The inlet pressure however is severely under-predicted because the pressure losses in the stenotic region are not accounted for correctly. The value of $\zeta=2$ corresponds to a Poisseuille profile while the DNS results of Varghese et al. (2007a, b) indicate that in the stenosis the profile becomes flatter, leading to much higher shear stress. For example, for steady inlet $Re=500$ , the wall shear stress varies very sharply across the stenotic region (refer to figure 18 of Varghese et al. (2007a)) and at the throat is $22$ times larger compared to upstream levels. Interestingly, the location of maximum shear amplification factor (equal to 30) is immediately prior of the throat. At $Re=1000$ , the maximum values increase by more than a factor of 40.

Based on this evidence from DNS results, we performed an additional computation in which $\zeta=44$ in the stenotic region and $2$ everywhere else. The large value of $\zeta$ corresponds to an average increase of the shear by $11.5$ times compared to Poisseuille profile. The results are shown in figure 11. It is clear that the peak value of pressure is now much better approximated but the secondary peak and the values during systole are slightly overpredicted. As $\zeta$ depends strongly on $Re$ , it is expected that more accurate predictions could have been obtained if a time- and space-dependent $\zeta$ is considered. However we have not pursued this further as the current results can be considered acceptable; for example the error at the peak is less than $3\%$ .

7 Conclusions

An 1D model to represent stenotic vessels has been proposed. Unlike models previously developed in the literature, the vessel is composed of series of tapered sections, instead of sections with constant area. At each section, the pressure and velocity can be expressed analytically in terms of Bessel functions. Using appropriate compatibility and boundary conditions between the sections, the full 1D pressure and velocity fields can be quickly obtained numerically in the frequency domain. The superposition of forward and reflecting waves results at resonant frequencies.

By representing the stenosis using only two segments, we can derive closed form analytical solutions that match very well the numerical solution at small, physiologically relevant, frequencies. The 1D model results were compared against 3D solutions that account for fluid-structure-interaction. Excellent matching was achieved between for the volume flow rate, but differences were noted on the pressure distribution. Using input from DNS simulations to approximate better the velocity profile around the stenotic region, significantly better results were obtained.

The proposed method can be used to approximate arbitrary profile shapes. It can also be extended to handle the whole cardiovascular tree. This will make sensitivity and uncertainty quantification studies much more affordable than they currently are.

Bibliography29

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Ahmed and Giddens (1983) Ahmed, S.A., Giddens, D.P., 1983. Velocity measurements in steady flow through axisymmetric stenoses at moderate Reynolds number. Journal of Biomechanics 16(7), 505––516.
2Alastruey et al. (2011) Alastruey, J., Khir, A.W., Matthys, K.S., Segers, P., Sherwin, S., Verdonck, P., Parker, K., Peiró, J., 2011. Pulse wave propagation in a model human arterial network: Assessment of 1-D visco-elastic simulations against in vitro measurements. Journal of Biomechanics 44(12), 2250–2258.
3Alastruey et al. (2012) Alastruey, J., Parker, K.H., Sherwin, S.J., 2012. Arterial pulse wave haemodynamics, in: 11th International Conference on Pressure Surges. Virtual Pi E Led t/a BHR Group, pp. 401–443.
4Berger and Jou (2000) Berger, S., Jou, L.D., 2000. Flows in stenotic vessels. Annual Review of Fluid Mechanics 32, 347–382.
5Cassanova and Giddens (1978) Cassanova, R.A., Giddens, D.P., 1978. Disorder distal to modeled stenoses in steady and pulsatile flow. Journal of Biomechanics 11, 441–453.
6Chen et al. (2013) Chen, P., Quarteroni, A., Rozza, G., 2013. Simulation-based uncertainty quantification of human arterial network hemodynamics. International Journal for Numerical Methods in Biomedical Engineering 29, 698–721.
7Choi et al. (2017) Choi, W., Park, S., Huh, H., S.J., L., 2017. Hemodynamic characteristics of flow around a deformable stenosis. Journal of Biomechanics 61, 216–223.
8Du et al. (2016) Du, T., Hu, D., Cai, D., 2016. A fast algorithm for the simulation of arterial pulse waves. Journal of Computational Physics 314, 450–464.