<figure><img src="image_1.png" /></figure>

HAL
open science

# Implementing Boundary Conditions in Simulations of Arterial Flows

Plamen Bokov, Patrice Flaud, Aziz Bensalah, Jose-Maria Fullana, Maurice Rossi

► To cite this version:

Plamen Bokov, Patrice Flaud, Aziz Bensalah, Jose-Maria Fullana, Maurice Rossi. Implementing Boundary Conditions in Simulations of Arterial Flows. Journal of Biomechanical Engineering, 2013, 135 (11), pp.111004. 10.1115/1.4025111. hal-01447079v2

HAL Id: hal-01447079

https://hal.science/hal-01447079v2

Submitted on 3 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Implementing Boundary Conditions in Simulations of Arterial Flows

Plamen Bokov, Patrice Flaud, A. Bensalah, Jose-Maria Fullana, Maurice Rossi

► To cite this version:

Plamen Bokov, Patrice Flaud, A. Bensalah, Jose-Maria Fullana, Maurice Rossi. Implementing Boundary Conditions in Simulations of Arterial Flows. 2013. hal-00800820

HAL Id: hal-00800820

https://hal.sorbonne-universite.fr/hal-00800820

Preprint submitted on 14 Mar 2013

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Implementing Boundary Conditions in Simulations of Arterial Flows

**P. Bokov and P. Flaud**

Laboratoire Matière et Systèmes Complexes
Université Paris Diderot
75251 Paris, France

**A. Bensalah**

CNRST BP 8027
10102 Rabat Maroc

**J.-M. Fullana and M. Rossi**

Institut Jean Le Rond D'Alembert
UPMC Univ Paris 06, CNRS UMR 7190
F-75005, Paris, France
Email: jose.fullana@upmc.fr

submitted 02/2013 - Journal of Biomechanical Engineering

Computational hemodynamic models of the cardiovas-
cular system are often limited to finite segments of the
system and therefore need well-controlled inlet and out-
let boundary conditions. Classical boundary conditions are
measured total pressure or flow rate imposed at inlet and
impedances of RLR, RLC or LR filters at outlet. We present
a new approach based on a unidirectional propagative hy-
pothesis (UPH) to model the inlet/outlet boundary condi-
tions on the axisymmetric Navier-Stokes equations. This
condition is equivalent to a non-reflecting boundary condi-
tion in a fluid-structure interaction model of an axisymmetric
artery. First we compare the UPH to the best impedance fil-
ter (RLC). Second, we apply this approach to a physiological
situation i.e. the presence of a stented segment into a coro-
nary artery. In that case a reflection index is defined which
quantifies the amount of pressure waves reflected upon the
singularity.

# 1 Introduction

As widely documented [ [1-8]], hemodynamical condi-
tions at the arterial wall play an important role in atheroscle-
rosis: a low mean shear stress, an increase of the shear oscil-
latory component, or an abnormal stagnation time are known
to be correlated to the location of atheromateous plaques. A
comprehensive study of such a phenomenon requires the ac-
curate knowledge of blood flow, i.e. velocity profiles and
pressure distributions inside the artery. Unfortunately, these
data remain difficult to measure experimentally, particularly
when studying atheromateous plaques: blood velocity mea-
surements using Doppler devices are difficult to assess, in

part because of the size and location of the region of inter-
est and, more importantly, because of the presence of distal
walls that deteriorates the signal to noise ratio due to wave
reflections. Another important point concerns the evaluation
of in vivo pressures, which cannot be directly obtained in an
atraumatic way. Because of these experimental limitations,
computational hemodynamics has been used as an alternative
tool for studying arterial flow. As large numerical resources
are required to simulate a complex 3D arterial tree, numerical
works are often limited to the investigation of a unique vascu-
lar segment [ 9–11]]. In all cases the main problem remains
the choice of inlet, outlet and arterial wall conditions. Differ-
ent approaches have been proposed to determine an optimal
choice for inlet and outlet conditions: for instance, [12] im-
posed inlet and outlet pressure to be the experimentally mea-
sured quantities. Resulting velocity profiles are nevertheless
sensitive to errors resulting from imprecise estimations of
pressure. Other authors used a RLC system (like Winkessel
models) or more complex systems with inductive elements
[ 13–16]]. These methods however failed to completely
eliminate inlet and outlet boundary reflections [ [17]]. This
latter point is crucial since spurious reflections generated by
outlet or inlet boundaries may flaw the investigation of true
reflections due to arterial stenosis or stents. In the follow-
ing, we present a numerical model of an axisymmetric flow
in a flexible tube and an approach based on a unidirectional
propagative hypothesis (UPH). This models the inlet/outlet
boundary conditions for the Navier-Stokes equations. We
compare the UPH with the condition obtained when assum-
ing a semi-infinite tube or a classical RCL system. Since
spurious reflections disappear in that approach we apply the
UPH to investigate the hydrodynamic role of a singularity;
more precisely we define a reflection index that quantifies the amount of pressure waves reflected upon a stent model.

## 2 Material and methods

We have tested the different boundary conditions on a Newtonian flow in a model of artery which is a flexible cylinder of length L and radius R(t). The temporal dependence on radius R(t) is due to the coupling between flow and wall. Thereafter, we present the numerical description of the model and its validation.

### 2.1 Description of the model

The numerical simulations were performed using the COMSOL software (COMSOL Inc.). This tool based on a finite element approach, numerically solves partial differential equations (PDE) with the corresponding boundary conditions. We use two different solvers one for the fluid domain and one for the solid one. COMSOL also handles fluid-structure problems flow by coupling both fluid and solid solvers. For the fluid part, COMSOL solves the Navier-Stokes equation and mass conservation over the computational domain. The fluid is assumed incompressible of constant density $\rho = 1060\,kg/m^3$ and Newtonian of dynamic viscosity $\mu = 5 \cdot 10^{-3}\,Pa \cdot s$. The blood viscosity for a high shear rate is around $4 \cdot 10^{-3}$ and can rise up to $40 \cdot 10^{-3}\,Pa \cdot s$ for low shear rates. We decided to use a moderate high viscosity value to model large or middle size arteries. The flow is assumed to remain axisymmetric; as a consequence, only radial $v_r(r,z,t)$ and longitudinal $v_z(r,z,t)$ velocities and pressure $p(z,t)$ appear into the mathematical model. The fluid equations thus read:

$$ 
\begin{aligned}
\frac{\partial v_r}{\partial t} + v_r \frac{\partial v_r}{\partial r} + v_z \frac{\partial v_z}{\partial z} &= \\
-\frac{1}{\rho} \frac{\partial p}{\partial r} + v \left( \frac{\partial^2 v_r}{\partial^2 r} + \frac{1}{r} \frac{\partial v_r}{\partial r} + \frac{\partial^2 v_z}{\partial^2 z} - \frac{v_r}{r^2} \right) \frac{\partial v_z}{\partial t} + v_r \frac{\partial v_z}{\partial r} + v_z \frac{\partial v_z}{\partial z} &= \\
-\frac{1}{\rho} \frac{\partial p}{\partial z} + v \left( \frac{\partial^2 v_z}{\partial^2 r} + \frac{1}{r} \frac{\partial v_z}{\partial r} + \frac{\partial^2 v_z}{\partial^2 z} \right) \frac{\partial v_r}{\partial r} + \frac{v_r}{r} \frac{\partial v_r}{\partial r} + \frac{v_z}{r} \frac{\partial v_z}{\partial r} &= 0
\end{aligned}
\quad (1) $$

where $v = \mu/\rho$ denotes the kinematic viscosity and the longitudinal variable $z$ runs between 0 and $L$, and the variable $r$ between 0 and $R(t)$. Furthermore, by symmetry, radial velocity and radial derivative of axial velocity are zero at $r = 0$ and a no slip condition is imposed on the arterial wall ($r = R(t)$).

In COMSOL, the inlet and outlet boundary conditions ($z = 0$ and $z = L$) are expressed as follows: the pressure is imposed at both ends of the computational domain. Depending on models (see below), it can be given directly as a given input ($P(z = 0/L, t)$), or can be obtained through the solution of an ordinary differential equation relating pressure and flow rate.

For the solid dynamics part, the arterial wall is assumed to move only radially which is generally well admitted in hemodynamics. The entrance and the outlet portion of the solid is free to move radially under the pressure forces. The flexible wall is assumed to be nearly incompressible and made of an isotropic linearly elastic solid. It is thus only defined by its Young Modulus $E$ and its Poisson's ratio $\nu_p$. The Euler equation of the solid structure is then:

$$ \frac{\partial \sigma_{ij}}{\partial x} = \rho \frac{d^2 b_i}{dt^2} \quad (2) $$

where $\sigma_{ij}$ denotes the component $(r, \theta, z)$ of the stress tensor and $b_i$ the component of deformation along the $i$-axis. For a linear elastic, isotropic material the stress tensor components are related to the strain tensor components $\varepsilon_{ij}$ by the relation:

$$ 
\begin{aligned}
\sigma_{ij} &= \frac{E}{1 + v_p} \left( \varepsilon_{ij} + \frac{v_p}{1 - 2 \cdot v_p} \varepsilon_{kk} \delta_{ij} \right), \\
\varepsilon_{ik} &= \frac{1}{2} (\partial_k b_i + \partial_i b_k + \partial_i b_l \partial_k b_l)
\end{aligned}
\quad (3) $$

where $\delta_{ij}$ is the Kronecker tensor. Using equations (2) and (3) COMSOL obtain a global equation for the displacement tensor. Equation (3) defines the two wall models used in our numerical computations. In the first one, the Young modulus $E$ is constant and, in the second one, $E = E(z)$ varies along the arterial axis ($z$ variable) to mimic the effect of the stent. This approach is valid: the Young modulus of the stent is known to be higher than the Young modulus of the arterial wall and, as shown in a previous work ([18]), the stent behaves as if it only induces a mechanical stiffening of the arterial wall. Practically, Young's modulus is increased for a normal value (artery) to a higher value within the length of the stent (20mm). In order to avoid discontinuous functions in numerical resolution the Young modulus smoothly changes near of the extremities of the stent.

To couple both solvers, COMSOL uses an ALE (Arbitrary Lagrangian-Eulerian) approach: the fluid-structure interaction is taken into account by imposing that radial forces computed through the fluid stress tensor are continuous across the wall at $r = R(t)$.

### 2.2 Computational domain

The computational mesh for the fluid domain is axisymmetric as stated above, therefore $(r, z = [0, R(t)] \times [0, L])$ is composed of rectangular elements: 500 elements in the axial direction and 12 elements in the radial one. We performed several sensitivity tests by simulating a pulsed flow in a long tube ($L = 0.5\,m$). For the axial direction, we studied the dependence of velocity field and arterial radius with the mesh density. The computation of the maximum radius over the simulation $Max_i[R(L/5, t)]$ at $z = L/5$ only varies of 0.04% when using 1000 instead of 500 axial elements. At the same position, the axial velocity varies by 1.7% and the wall
shear stress (WSS) by 1.9%. On the contrary, 250 mesh elements seem to be insufficient. Therefore for all numerical simulations presented here the mesh density in the axial direction was set to 500/50 = 10 elements per cm. This radial mesh is not uniform: step size is growing exponentially as $\exp(0.03r)$ from the wall to the center of the tube. This allows meshing correctly the oscillating boundary layer; its typical size is proportional to the inverse of the square of the Womersley number $\alpha$ defined as $R(\omega/v)^{1/2}$. Here $\omega$ denotes the typical flow pulsation.

We have checked that 12 elements in the radial direction were sufficient enough for a good evaluation of the WSS. For the elastic wall, only 3 mesh elements were used in the radial direction. We tested the sensitivity of the solution and determined that the relative variation of the solution between 2 elements and 5 elements in the tube thickness was less than 0.3% for the maximal fluid velocity at $z = L/5$. The axial mesh in the solid is identical to the fluid domain one. Consequently, the sensitivity with respect to the axial mesh in the solid has been performed together with the fluid domain (see above).

## 2.3 Validation of the model

In order to validate our numerical model, we compare the numerical solution with an analytical solution in a semi-infinite medium. Note that this analytical solution is only pertinent in a linear regime where deformations amplitudes are small compared to physiological ones.

### Analytical solution in a "semi-infinite" tube

It is possible to obtain an analytical solution for a pressure wave generated by an imposed periodic pressure signal $p(z = 0, t)$ which propagates in a viscous fluid contained in a semi-infinite elastic tube $z > 0$. Indeed let us Fourier decompose the pressure at $z = 0$

$$P(t, z=0) = \sum_{n=0}^{N} p_n e^{i\omega_n t}.$$

The general solution then reads

$$P(t, z) = \sum_{n=0}^{N} p_n e^{i(\omega_n t - k_n z)}$$

where the complex wave numbers $k_n(\omega_n)$ is obtained by a dispersion relation for a particular propagative solutions of the form $e^{i(\omega t - k(\omega)z)}$. Within some approximations (wave velocity larger than the longitudinal fluid velocity, wavelength larger than the mean tube radius $R$), the existence of this particular solution implies that complex wave number $k$ and the real pulsation $\omega$ are related via the dispersion equation

$$c^2(\omega) = \frac{\omega^2}{k^2} = c_0^2 \left[ 1 - F(i^{3/2}\alpha) \right], \quad \alpha = R \left[ \frac{\omega}{v} \right]^{1/2} \quad (4)$$

where $F(z) = 2J_1(z)/(zJ_0(z))$ is the Womersley function [19], $J_0$ and $J_1$ are Bessel functions of the 1st kind, and $c_0$ is the phase velocity of the radial mode called Young's mode equal to $\frac{A}{\rho} \frac{dP}{dA}$. For thin walled tube hypothesis and linear elasticity, its value is known via the Moens and Korteweg relation [20]

$$c_0^2 = \frac{A_0}{\rho C_p} \quad (5)$$

where $A_0$ stands for the undeformed tube cross sectional area and $C_p$ for the elastic tube compliance. This latter quantity is assumed constant (Hookean model) which is the correct approximation as long as the pressure-induced deformations are negligible compared to the tube radius.

### Numerical simulation in a "semi-infinite" tube

The above analytical expression is used to check the numerical code for a long tube with a known inlet pressure. The pressure imposed at inlet $P_{imp}(t)$ mimics a systolic pulse of duration $T_1 = 0.4$ s and of period $T = 1$ s, then

$$P_{imp}(t) = \begin{cases} \Delta P (1 - \cos(2\pi t/T_1)) & t \in [0, T_1] \\ 0.0 & t \in [T_1, T] \end{cases} \quad (6)$$

where $\Delta P$ is the maximal pressure pulse.

<table><caption>Table 1. Mechanical properties and numerical inputs used for the "semi-infinite" simulations.</caption><tr><td>h/D<sub>ext</sub></td><td>0.083</td></tr><tr><td>Tube length [m]</td><td>4.92</td></tr><tr><td>ΔP/E</td><td>5.10<sup>-4</sup></td></tr><tr><td>Womersley number α</td><td>8.2</td></tr><tr><td>Dynamic viscosity (blood) [Pa.s]</td><td>5.10<sup>-3</sup></td></tr><tr><td>Pulse duration T<sub>1</sub> [s]</td><td>0.4</td></tr><tr><td>Computed attenuation coefficient [m<sup>-1</sup>]</td><td>0.20</td></tr></table>

Mechanical and geometrical properties as well as numerical inputs for COMSOL are given in Table 1 where $h$ denotes wall thickness, $D_{ext}$ the external diameter.

These values are chosen close to physiological values of human coronary arteries since our interest is focused towards studying pressure pulse propagation in coronary arteries and the effect of stent placement in such vessels. Physiological $h/D_{ext}$ for coronary arteries have been extensively studied: for instance, [21] reported values between 0.08 and 0.44 depending on the age of the subject and the artery (left coronary artery vs right coronary artery).

The simulation is called here "semi infinite" because the tube is long enough ($L = 4.92m$) to avoid wave reflections at
the outlet over the simulation time. In Figure 1, comparisons are performed for $z \in [0, 1.6 m]$ and reflected waves have no time to come-back to measurement points and only one wave is observed, the incident wave.

Comparisons (Figure 1) confirms the validity of our numerical code. The pulse attenuation coefficient was estimated by fitting the decrease of peak-to-peak pressure by an exponential function of the distance. The axial velocity of the fluid is of the order of $3.10^{-2}$ m/s and the wave group velocity (peak to peak) computed from the numerical simulation was $5.7 \pm 0.1$ m/s (the Moens and Korteweg velocity (equation (5)) is 6.14 m/s which corresponds to the high frequency phase velocity). In addition, linear hypotheses are fulfilled since non-linear terms are around 0.5% of the linear term $\frac{\partial u}{\partial t}$.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Non dimensional pressure (P/E) as a function of time. Each curve corresponds to different positions z (0, 0.4m, 0.8m, 1.2m, and 1.6m). Numerical simulations (red dash-dot line) are compared to analytical solutions (black dashed line). The peak Reynolds number evaluated at entry, is equal to 23.</figcaption></figure>

## 3 Boundary Conditions

In this section, we propose a new boundary condition to allow unidirectional waves cross boundaries without reflections, the approach is named “unidirectional propagation hypothesis” (UPH). We also build another boundary condition based on a RCL model. Both approaches will be compared in Section 4.

### 3.1 Outlet boundary condition: unidirectional propagation hypothesis (UPH)

Let us consider the axial-component of the linearized Navier-Stokes equation in which one assumes axial derivatives to be negligible with respect to radial derivatives (such hypotheses are justified for blood flow since artery radius is smaller enough than characteristic wavelength). This leads to the simplified equation

$$ \frac{\partial v_z}{\partial t} = -\frac{1}{\rho} \frac{\partial P}{\partial z} + \mu \left[ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial v_z}{\partial r} \right) \right] \quad (7) $$

where P is the transmural pressure. When integrated over a section at constant z, this equation becomes:

$$ \frac{dQ}{dt} = -\pi R^2 \frac{1}{\rho} \frac{\partial P}{\partial z} + 2\pi R \mu \frac{\partial v_z}{\partial r} \Big|_{r=R} \quad (8) $$

Let us consider a monochromatic wave of pulsation $\omega$ propagating in the $z > 0$ direction. Within Prandtl approximation, the associated transmural pressure is independent of the radial coordinate r and tis a plane progressive wave damped in the propagation direction:

$$ P(z,t) = P_0 e^{\gamma z} e^{i\omega(t-z/c)} \quad (9) $$

where $c = c(\omega)$ is the phase velocity and $\gamma = \gamma(\omega)$ a damping coefficient. For this particular solution propagating towards positive z, it is possible to transform the z derivative into a time derivative

$$ \frac{\partial P}{\partial z} = -\gamma P - \frac{1}{c} \frac{\partial P}{\partial t} \quad (10) $$

Let us now assume that this relation can be used for any unidirectional pressure wave as long as c and $\gamma$ only a slightly depend on frequency $\omega$. Using equations (8) and (10), a relationship is obtained for outlet pressure $P_{out} = P(z_{out}, t)$:

$$ -\gamma P_{out} + \frac{1}{c} \frac{dP_{out}}{dt} - \frac{\rho}{\pi R^2} \frac{dQ_{out}}{dt} + \frac{2\mu}{R} \frac{\partial v_z}{\partial r} \Big|_{r=R, z=z_{out}} = 0 \quad (11) $$

In our computations, this equation is implemented as the outlet condition. Indeed pressure $P_{out}(t)$ can be obtained by solving this equation, $Q_{out}$ being numerically computed. The following procedure is called the unidirectional propagation hypothesis (UPH).

The values of c and $\gamma$ used in equation (11) are determined through simulations as follows. First a computation is performed using outlet condition (equation (11)) with constant c given by the Moens-Korteweg relation (5) and $\gamma$ by the attenuation in the case of Womersley number given in Table 1. Moreover, we use for this first simulation a tube long enough (1 m in this case) to ensure the separation of backward and forward waves. By following the position and value of the pulse peak, velocity c and attenuation $\gamma$ are measured. Once these factors computed, a second simulation is performed with the new values of c and $\gamma$ and the true finite tube length. If reflected waves are observed, one repeats the same procedure and computes again the values of c and
$\gamma$ until the reflective waves disappear. Convergence is generally reached in three iterations. The evaluation of parameter $c$ is important since a wrong value of $c$ would lead to an unsatisfactory non-reflective condition. If one wants to modulate a reflected wave by prescribing a given amplitude index $a_{in}$ which runs between 0 (no reflexion) and 1 (total reflexion), it is possible to use the modified equation

$$a_{in} \left( -\gamma P_{out} + \frac{1}{c} \frac{dP_{out}}{dt} \right) - \frac{\rho}{\pi R^2} \frac{dQ_{out}}{dt} + \frac{2\mu}{R} \frac{\partial v_z}{\partial r} \bigg|_{r=R,z=z_{out}} = 0$$

It is possible to introduce similar boundary conditions at inlet to minimize reflective waves generated at the entrance. It is however necessary to separate the total pressure as a sum of a wave $P_i(z,t)$ propagating in the increasing $z$ direction and a wave $P_r(z,t)$ propagating in the decreasing $z$ direction:

$$P(z,t) = P_i(z,t) + P_r(z,t) \quad (12)$$

with $P_i(z=0,t) = P_{imp}(t)$. One can then write equation (8) for pressure $P$ and flow rate at inlet ($z=0$). The reflected wave $P_r$ and the imposed pressure function $P_{imp}$ at that point can be added to get the total pressure $P_{in} = P(z=0,t)$. Assuming that $P_i(z,t)$ is given by a spatio-temporal function similar to equation (9) and that $P_r(z,t)$ is given by

$$P_r(z,t) = P_{r0} e^{\gamma z} e^{i\omega(t-z/c)}$$

then the $z$-derivatives of $P(z,t)$ at the entrance read:

$$\begin{aligned} \frac{\partial P_{in}(z=0,t)}{\partial z} &= \frac{\partial P_i(z=0,t)}{\partial z} + \frac{\partial P_r(z=0,t)}{\partial z} \\ &= -\gamma P_{imp} - \frac{1}{c} \frac{dP_{imp}}{dt} + \left[ \gamma P_r + \frac{1}{c} \frac{dP_r}{dt} \right]_{z=0} \end{aligned} \quad (13)$$

Using equation (12), one can eliminate $P_r$ from equation (13) and deduce a boundary condition at the entrance:

$$\left( \gamma P_{in} + \frac{1}{c} \frac{\partial P_{in}}{\partial t} \right) - \gamma P_{imp} - \frac{1}{c} \frac{\partial P_{imp}}{\partial t} + \frac{\rho}{\pi R^2} \frac{dQ_{out}}{dt} + \frac{2\mu}{R} \frac{\partial v_z}{\partial r} \bigg|_{r=R} = 0$$

Equations (11) and (14) are used as boundary conditions: they allow to model pressure pulses with minimal reflections at both tube extremities. This model has one major assumption: the non-dispersive hypothesis concerning $c$ and $\gamma$.

A priori the viscous term in equation (14) must be evaluated at $z=0$. Nevertheless at the inlet, the flow is not fully developed and the wall shear rate computed at $z=0$ would lead to a wrong estimation of this viscous term. Under the long wavelength hypothesis it is more convenient to evaluate this term at a position sufficiently away from the inlet, larger than the inlet length and negligible as compared to the wavelength. The extent along $z$ of this region scales as $R$,

the radius of the tube. The effect of an inaccurate evaluation of the near wall shear rate can be observed by considering a wave propagating towards the inlet as shown in Figure 2a. If one evaluates the shear rate of equation (14) in a point immediately near the entrance (point C in Figure 2a) one can observe reflections (Figure 2c). Furthermore when the evaluation is done in a point distant from the inlet (point B located at a distance of $7R$ from the inlet), no reflective wave appears (Figure 2b).

<figure><img src="image_4.png" /><figcaption>Fig. 2. (a) Pulse wave, propagating towards negative $z$: axial velocity is color coded (max in red, deep blue is zero) and the contour is the radial displacement for three times. (b) Pressure as a function of time at two different positions, $z=0.4m$ and $z=0.1m$. The shear rate of equation (23) is evaluated at $z=1cm$, i.e. approximately $7R$, (c). Idem Fig 2 (b) but the shear rate is now evaluated at the tube entrance.</figcaption></figure>

## 3.2 Another outlet boundary condition : RLC as a model of terminal impedance

We present in this section a boundary condition based on a RCL filter model. If one assumes low amplitude, the outlet pressure is related to the flow rate by the characteristic impedance $Z_c(\omega)$

$$Z_c(\omega) = \frac{\rho c_0}{A_0} \left[ 1 - F(\dot{t}^{3/2}\alpha) \right]^{-1/2} \quad (15)$$

This relation however cannot be used with a time advancing computation scheme since, it requires to work in frequency spectral space and thus to know, a priori, the whole time signal. In the general nonlinear regime, we replaced this transformation by a differential equation system approximating the function $Z_c$ in a given frequency range, in particular the physiological one, between 1 and 20 Hz. We propose a more complex filter than a standard RLC, where the elements are estimated by fitting the filter impedance as closely as possible to the characteristic impedance. When linear hypothesis are fulfilled, such a condition is strictly equivalent
to the semi-infinite tube. This particular filter is given by a resistance $R_1$ and an inductive element $L$ in series with a capacity $C$ in parallel with an element $Z_2$. This leads to the following complex impedance

$$Z = R_1 + i\omega L + \frac{Z_2}{1 + i\omega C Z_2}$$

In order to obtain the best fit of equation (15), we assumed that element $Z_2$ is complex. Once the best values of coefficient are found ($C, R_1, L, Z_2$) it is possible to connect such an impedance to a system of differential equations calculating the pressure/flow relationship at the outlet boundary. In Fourier frequency-space, the expression for $Z$ is decomposed into two terms $Z = Z_1 + 1/Y_2$, so that the total complex pressure $\bar{P}$ at the outlet can be split into to terms such $\bar{P} = \bar{P}_1 + \bar{P}_2$

$$\bar{P}_1 = Z_1 Q = (R_1 + i\omega L) Q \\ \bar{P}_2 = \frac{1}{Y_2} Q = \frac{1}{\left(\frac{1}{Z_2} + i\omega C\right)} Q$$

where $Q$ is the real flow rate. The above algebraic system is equivalent to a system of linear differential equations which will be used as an implicit outlet boundary condition:

$$\bar{P}_1 - R_1 Q - L \frac{dQ}{dt} = 0 \\ \frac{\bar{P}_2}{Z_2} + C \frac{d\bar{P}_2}{dt} - Q = 0.$$

Finally the pressure $P_{out}$ is obtained by the following:

$$P_{out} = \mathcal{R}_e(\bar{P}_1 + \bar{P}_2).$$

In reference [18], we have shown that this approach is better than the classical RLC filter method to define non reflective boundary conditions.

## 4 Results and Discussion

We present numerical test on the application of the UPH approach. Let us denote by $S_{INF}$ the value of the numerical solution of a pulse wave propagating in a semi-infinite domain (see section 2) but reduced to the spatial sector $[0, 0.5m]$. At $z = 0.5m$, the numerical data hence present no reflexion and no reflected waves. Therefore this numerical data is the target for our comparisons. $S_{UPH}$ denotes the solution obtained numerically within domain $[0, L = 0.5m]$ using the UPH boundary condition at $z = 0.5$ and $S_{RLC}$ the associated solution obtained numerically in the domain $[0, L = 0.5m]$ using the RLC outlet boundary condition. We simulate these various approximations to identify which boundary condition is preferable. The mechanical characteristics and simulation parameters are given in Table 2.

## 4.1 Comparison between UPH versus RLC.

<table><caption>Table 2. Mechanical properties and numerical input used in UPH simulations.</caption><tr><td>Internal radius [mm]</td><td>1.5</td></tr><tr><td>Tube length [m]</td><td>0.5</td></tr><tr><td>h/D<sub>ext</sub></td><td>0.083</td></tr><tr><td>Pulse period [s]</td><td>0.075</td></tr><tr><td>ΔP/E</td><td>5.10<sup>-4</sup> or 5.10<sup>-3</sup></td></tr><tr><td>Moens-Korteweg velocity [m.s<sup>-1</sup>]</td><td>6.32</td></tr><tr><td>Womersley number (α)</td><td>6.14</td></tr><tr><td>Computed group velocity [m.s<sup>-1</sup>]</td><td>5.87 ± 0.05</td></tr></table>

We use the imposed pressure signal $P_{imp}$ of equation (6) with two values of adimensional pulse pressure amplitude $\Delta P/E$ ($5.10^{-4}$ and $5.10^{-3}$) which represents approximately 1% and 10% of physiological values. The last value corresponds to the limit of approximations in the analytical model as $v_{max}/c$ becomes close to 6%.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Adimensional pressure ($P/E$) as a function of time at several axial positions ($z=0$ m (entry point); 0.1m; 0.3m; 0.5m (exit point)). Red dash-dot lines correspond to $S_{RLC}$, black dashed lines to $S_{INF}$. (a) $\Delta P/E$ of $5.10^{-4}$; (b) $\Delta P/E$ of $5.10^{-4}$.</figcaption></figure>

For $\Delta P/E = 5.10^{-4}$, the two solutions $S_{RLC}$ to $S_{INF}$ agree (Figure 3a), the relative difference being less than 3.5%. For $\Delta P/E = 5.10^{-3}$ (Figure 3b), the difference at the tube outlet grows to 7% due to reflections.

For $\Delta P/E = 5.10^{-4}$, the solution $S_{UPH}$ and $S_{INF}$ are in good agreement (Figure 4a). The wave velocity and attenuation coefficient are also well assessed. For $\Delta P/E = 5.10^{-3}$ (Figure 4b), the difference between the two solutions at the outlet boundary is around 0.25% (evaluated at pressure peak). This result has to be compared to the $S_{RLC}$ condition (7%).

Figure 2a illustrates the rôle of the inlet boundary condition (14): in that instance, a pressure wave propagates towards the inlet and exists without reflections. However, when both backward and forward waves interact at $z = 0$, this inlet condition has not as efficient.
The UPH boundary conditions at inlet and outlet is able to eliminate spurious reflective waves. As illustrated in Figure 4, it provides for a short tube in a 2D-axisymmetric configuration, a pressure propagative wave which does not differ from the pressure computed in “infinite” domain by more than 0.25%.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Non-dimensional pressure ($P/E$) as function of the time at several axial positions (z=0 m (entry point); 0.1 m; 0.3 m; 0.5 m (exit point)). Red dash-dot lines correspond to $S_{UPH}$ with UPH outlet conditions (equation 18), black dashed lines to $S_{INF}$ (a) $S_{UPH}$ and $S_{INF}$ for $\Delta P/E = 5.10^{-4}$ (b) $S_{UPH}$ and $S_{INF}$ for $\Delta P/E = 5.10^{-3}$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Pressure pulse using UPH for a high amplitude ($10^4$ Pa) in a 50 cm length tube with Young modulus $E = 106$ Pa and viscosity $\mu = 7$ mPa.s. Wall displacements at successive times (0.05, 0.1, 0.15, 0.2, 0.25 and 0.3 sec) are shown: no apparent reflexion occurs. The wave moves from the bottom to the top. The color bar is scaled in $10^{-2}$ mm, and the bottom axis in mm.</figcaption></figure>

## 4.2 Application of UPH in high amplitude pressure pulse

The above result (Figure 4 a) is still valid when $\Delta P/E$ is increased by an order of magnitude (Figure 4 b) which induces, for a given Youngs modulus, wall displacements about 10 times larger. For even larger $\Delta P/E$ values, which are close to characteristic physiological values, this conclusion remains valid. More precisely, we tested the efficiency of the UPH for a $10^4$ Pa amplitude in a 0.5 m tube length. This amplitude is of the same order of magnitude of the systolic pulse, i.e. the difference between the systolic and diastolic pressures encountered in human circulation during the systole. The pulse period (100 ms) is comparable to the systolic rise time. Figure 5 presents snapshots of the wall displacement for times 0.05, 0.1, 0.15, 0.2, 0.25 and 0.3 sec. The wave moves from the bottom to the top and exits without noticeable reflections. In a more quantitative way, Figure 6 indicates that the results obtained by the UPH method are similar to those of a long tube simulation $S_{INF}$.

<figure><img src="image_8.png" /><figcaption>Fig. 6. Pressure pulse using UPH for a high amplitude ($10^4$ Pa) : the position of the pressure maximum as a function of time. Diamonds ($S_{INF}$), squares ($S_{UPH}$).</figcaption></figure>

## 4.3 Application of UPH on a stent singularity

Diseases of the arterial system (aneurism, stenosis...) and surgical treatments (stent-grafts, vascular substitute, stent) induce modifications of mechanical and geometrical properties of the vascular wall. For instance the consequences of a stent in an artery have been widely studied (struts protrusion [22]), distribution of mechanical stresses on the artery [23, 24]). In particular it is known that a typical metal stent is around 5 times less compliant than the arterial wall [25]. One practical application of the UPH approach is to describe quantitatively the hydrodynamic effect of such an inhomogeneity in a short vessel segment. In this case it is necessary to separate the induced reflected wave

from the total pulse wave (a stent could be responsible for a pressure pulse reflection around 6%).

Numerous separation methods have been presented in the literature [26, 27]: linear methods, which give fairly good results, however fail when considering high amplitude waves [18] and nonlinear methods used until now generally required strong hypothesis for the evaluation of viscous term in the Navier-Stokes equation. The proposed approach evaluates the impact of an inhomogeneity by (a) performing two numerical simulations : without and with inhomogeneity. and (b) by subtracting both solutions, upstream of the singularity. Since reflected waves are only issued from the singularity, the difference between the two solution signals (flow rate, pressure, shear stress...) gives the reflected wave.

Figure 7 presents successive wall displacement profile $\Delta R(z)$ in the tube measured each 1 m.sec. Without any singularity (fig.7b), the envelop of the maximum wall displacements (yellow line) shows the expected attenuation of pressure maxima while the wave is propagating. When a stent is present, this curve exhibits an increase just before the singularity, associated to a backward reflected wave. The characteristics of the numerical simulation are given in Table 2
<figure><img src="image_9.png" /><figcaption>Fig. 7. Radial displacement waves measured every $1e^{-3}$ s using UPH: Fig. (a) when the tube presents a singularity a hundred times stiffer than the tube wall situated between $z = 0.31$ m and $z = 0.33$ m; Fig. (b) without singularity. A yellow continuous line is superimposed, joining the maximum of the displacement wave at each time.</figcaption></figure>

with a pressure signal defined before, and $\Delta P/E$ of $5 \times 10^{-3}$. We have validated this method of separation by comparing its results within the linear domain with another separation method described in [27]. A good agreement was observed (results are not shown; differences are less than 0.5%).

A reflection index is introduced to quantify the effect of the singularity: it is based on the ratio between the backward axial velocity amplitude $W_{cl}$ and the forward amplitude $W_{cl}^+$ of the axial velocity $W_{cl} = w(r=0)$. Note that we also use an index based but on the backward tube deformation $\Delta R^-$ and forward amplitude $\Delta R^+$. The results are identical $W_{cl}^-; \Delta R^-$ and the forward amplitude $W_{cl}^+; \Delta R^+$ of the axial velocity $W_{cl} = w(r=0)$ and a second one is tube deformation induced by waves:

$$r = \frac{W_{cl}^-}{W_{cl}^+} \sim \frac{\Delta R^-}{\Delta R^+}$$

Figure 8 displays the numerical evaluation of this reflection index as a function of the relative Young's modulus of the stent.

The difference observed between a separation based on a linear theory of pulse reflection and the numerical separation method can be explained by the fact that a linear theory does not take into account the hydrodynamic entrance length in the stent. The reflection index is an indirect evaluation of the compliance of the stent. Such an evaluation could be interesting in monitoring the degradation of biodegradable stents. Indeed it is still unknown how long a biodegradable scaffold effectively assumes a supporting role of the vascular wall [28].

## 5 Conclusion

We introduced a method that defines a non-reflective boundary condition for pressure waves into a flexible vessel. This method takes also into account the fluid-structure interaction and is based on Navier-Stokes equations using inlet/outlet boundary conditions with an unidirectional propagative hypothesis (UPH). This approach allows to eliminate

<figure><img src="image_10.png" /><figcaption>Fig. 8. Reflection index as a function of the relative Young's modulus of the stent ($E_r = E/E_0$) where $E_0$ is the Young modulus of a normal artery. The results using our numerical separation method are given by red squares. The results of a linear theory of reflection (see Annex I) are represented by black squares.</figcaption></figure>

spurious reflections which appear in classical techniques like RLC and LR circuits. Though it is performed here with an axisymmetric geometry, this method can be easily generalized to 3D simulations. A reflection index is given by the implementation of a method of wave separation where the component of the reflected wave is computed from the values of total pressures or wall displacements with and without the obstacle (the stent). We have applied this approach to a physiological situation i.e. the presence of a stent into a coronary artery.

## References

[1] Nerem, R., and Cornhill, J., 1980. "The role of fluid mechanics in atherogenesis". *Journal of Biomechanical Engineering* (102), pp. 181-189.

[2] Zarins, C., Giddens, D., Bharadvaj, B., Sottiurai, V., Mabon, R., and Glagov, S., 1983. "Carotid bifurcation atherosclerosis. quantitative correlation of plaque localization with flow velocity profiles and wall shear stress". *Circulation Research* (53), pp. 502-514.

[3] Ku, D., Giddens, D., Zarins, C., and Glagov, S., 1985. "Pulsatile flow and atherosclerosis in the human carotid bifurcation. positive correlation between plaque location and low oscillating shear stress. atherosclerosis". *Thrombosis and Vascular Biology* (5), pp. 293-302.

[4] Giddens, D., Zarins, C., and Glagov, S., 1993. "The role of fluid mechanics in the localization and detection of atherosclerosis". *Journal of Biomechanical Engineering* (115), pp. 588-594.

[5] Tropea, B., Schwarzacher, S., Chang, A., Asvar, C., Huie, P., Sibley, R., and Zarins, C., 2000. "Reduction of aortic wall motion inhibits hypertension-mediated experimental atherosclerosis". *Arteriosclerosis, Thrombosis and Vascular Biology* (20), pp. 2127-2133.
[6] Meyerson, S., Skelly, C., Curi, M., Shakur, U., Vosicky, J., Glagov, S., and Schwartz, L., 2001. "The effects of extremely low shear stress on cellular proliferation and neointimal thickening in the failing bypass graft". *Journal of Vascular Surgery* (34), pp. 90–97.

[7] Cheng, C., Tempel, D., van Haperen, van der Baan, A., Grosveld, F., Daemen, M., Krams, R., and de Crom, R., 2006. "Atherosclerotic lesion size and vulnerability are determined by patterns of fluid shear stress.". *Circulation* (113), pp. 2744–2753.

[8] Coskun, A., Chen, C., Stone, P., and Feldman, C., 2006. "Computational fluid dynamics tools can be used to predict the progression of coronary artery disease". *Physica A* (362), pp. 182–190.

[9] Botnar, R., Rappitsch, G., Scheidegger, M., Liepsch, D., Perktold, K., and Boesiger, P., 2000. "Hemodynamics in the carotid artery bifurcation: a comparison between numerical simulations and in vitro mri measurements.". *Journal of Biomechanics* (33), pp. 137–144.

[10] Tang, D., Yang, C., Walker, H., Kobayashi, S., and Ku, D., 2002. "Simulating cyclic artery compression using a 3d unsteady model with fluid-structure interactions". *Computers & Structures* (80), pp. 1651–1665.

[11] Cebral, J., Castro, M., Appanaboyina, S., Putman, C., Milan, D., and Frangi, A., 2005. "Efficient pipeline for image-based patients specific analysis of cerebral aneurysm hemodynamics: technique and sensitivity". *IEEE Transactions on Medical Imaging* (24), pp. 457–467.

[12] Li, M., Beech-Brandt, J., John, L., Hoskins, P., and Eason, W., 2007. "Numerical analysis of pulsatile blood flow and vessel wall mechanics in different degrees of stenoses". *Journal of Biomechanics* (40), pp. 3715–3724.

[13] Vignon, I., and Taylor, C., 2004. "Outflow boundary conditions for one-dimensional finite element modeling of blood flow and pressure wave in arteries". *Wave Motion* (39), pp. 361–374.

[14] Vignon-Clementel, I., Figueroa, C., Jansen, K., and Taylor, C., 2006. "Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries". *Computer Methods in Applied Mechanics and Engineering* (195), pp. 3776–3796.

[15] Alastruey, J., Parker, K., Peiro, J., Byrd, S., and Sherwin, S., 2007. "Modelling the circle of williw to assess the effects of anatomical variations and occlusions on cerebral flows". *Journal of Biomechanics* (40), pp. 1794–1805.

[16] Kim, H., Vignon-Clementel, I., Coogan, J., Figueroa, C., Jansen, K., and Taylor, C., 2010. "Patient-specific modeling of blood flow and pressure in human coronary arteries". *Annals of Biomedical Engineering*, **10**(38), pp. 3195–3209.

[17] Willemet, M., Lacroix, V., and Marchandise, E., 2011. "Inlet boundary conditions for blood flow simulations in truncated arterial networks". *Journal of Biomechanics* (44), pp. 897–903.

[18] Bokov, P., 2011. "Description expérimentale et numérique de l'interaction entre un stent biodégradable et la paroi artérielle". PhD thesis, Université Paris Diderot.

[19] Bensalah, A., 1998. "Modèles d’écoulements artériels et quantification de la vélocimétrie ultrasonore a effet doppler pulsé en vue d’une étude atraumatique du système artériel". Master’s thesis, Université Denis Diderot.

[20] McDonald, 1960. *Blood Flow in Arteries*. Williams and Wilkings Company,.

[21] Ozolanta, I., Tetere, G., Purinya, B., and Kasyanov, V., 1998. "Changes in the mechanical properties, biochemical contents and wall structure of the human coronary arteries with age and sex". *Medical Engineering and Physics* (20), pp. 523–533.

[22] Seo, T., Schachter, L., and Barakat, A., 2005. "Computational study of fluid mechanical disturbance induced by endovascular stents". *Annals of Biomedical Engineering*, **33**(4), pp. 444–456.

[23] LaDisa, J.F, J., Olson, L., Guler, I., Hettrick, D., Kersten, J., Warltier, D., and Pagel, P., 2005. "Circumferential vascular deformation after stent implantation alters wall shear stress evaluated with time-dependent 3d computational fluid dynamics models". *Journal of Applied Physiology*:(98), pp. 947–957.

[24] Gijsen, F., Migliavacca, F., Schievano, S., Socci, L., Petrini, L., Thury, A., Wentzel, J., van der Steen, A., Serruys, P., and Dubini, G., 2008. "Simulation of stent deployment in a realistic human coronary artery". *BioMedical Engineering OnLine* (7), p. 23.

[25] Nicoud, F., Vernhet, H., and Dauzat, M., 2005. "A numerical assessment of wall shear stress changes after endovascular stenting". *Journal of Biomechanics* (38), pp. 2019–2027.

[26] Stergiopulos, N., Tardy, Y., and Meister, J., 1993. "Nonlinear separation of forward and backward running waves in elastic conduits". *Journal of Biomechanics* (26), pp. 201–209.

[27] Rogova, I., 1998. "Propagation d'ondes en hémodynamique artérielle: application à l'évaluation indirecte des paramètres physiopathologiques". PhD thesis, Université Paris Diderot. Paris. 1998.

[28] Wong, P., 2011. *Bioabsorbable scaffold: animal and human observations: Take-home message*.