A Novel Estimation Approach of Pressure Gradient and Haemodynamic Stresses as Indicators of
Pathological Aortic Flow Using Subvoxel Modelling
Pascal Corso, George Giannakopoulos, Utku G¨ulan, Christos Emmanouil Frouzakis, and Markus Holzner
Abstract—Objective: The flow downstream from aortic stenoses
1
is characterised by the onset of shear-induced turbulence that
2
leads to irreversible pressure losses. These extra losses represent
3
an increased resistance that impacts cardiac efficiency. A novel
4
approach is suggested in this study to accurately evaluate
5
the pressure gradient profile along the aorta centreline using
6
modelling of haemodynamic stress at scales that are smaller
7
than the typical resolution achieved in experiments.Methods: We
8
use benchmark data obtained from direct numerical simulation
9
(DNS) along with results from in silico and in vitro three-
10
dimensional particle tracking velocimetry (3D-PTV) at three
11
voxel sizes, namely 750 µm, 1 mm and 1.5 mm. A differential
12
equation is derived for the pressure gradient, and the subvoxel-
13
scale (SVS) stresses are closed using the Smagorinsky and a
14
new refined model. Model constants are optimised using DNS
15
and in silico PTV data and validated based on pulsatile in
16
vitro3D-PTV data and pressure catheter measurements.Results:
17
The Smagorinsky-based model was found to be more accurate
18
for SVS stress estimation but also more sensitive to errors
19
especially at lower resolution, whereas the new model was found
20
to more accurately estimate the projected pressure gradient even
21
for larger voxel size of 1.5mm albeit at the cost of increased
22
sensitivity at this voxel size. A comparison with other methods
23
in the literature shows that the new approach applied to in
24
vitro PTV measurements estimates the irreversible pressure
25
drop by decreasing the errors by at least 20%. Conclusion:
26
Our novel approach based on the modelling of subvoxel stress
27
offers a validated and more accurate way to estimate pressure
28
gradient, irreversible pressure loss and SVS stress. Significance:
29
We anticipate that the approach may potentially be applied to
30
image-basedin vivo,in vitro4D flow data orin silicodata with
31
limited spatial resolution to assess pressure loss and SVS stresses
32
in disturbed aortic blood flow.
33
Index Terms—aortic stenosis, turbulence, pressure drop, SVS
34
stress, DNS, 3D-PTV.
35
I. INTRODUCTION
36
A
ORTIC STENOSIS (AS) is a narrowing of the aortic37
valve area due to a progressively dysfunctional valve and
38
is the most common valvular disease in Europe and America
39
[27]. Increased resistance to the circulation of blood through a
40
P. Corso is with the Institute for Environmental Engineering, Swiss Federal Institute of Technology Zurich, Zurich, Switzerland (e-mail:
corso@ifu.baug.ethz.ch)
G. Giannakopoulos and Ch. E. Frouzakis are with the Aerothermochemistry and Combustion Systems Laboratory, Swiss Federal Institute of Technology Zurich, Zurich, Switzerland.
U. G¨ulan is with Hi-D Imaging, Zurich, Switzerland
M. Holzner is with the Swiss Federal Institute for Forest, Snow and Landscape Research WSL and the Swiss Federal Institute of Aquatic Science and Technology Eawag, Zurich, Switzerland.
high transvalvular pressure difference results in an increased 1
load on the left ventricle that can lead to hypertrophy and 2 dysfunction thereof. The prevalence of AS increases with age 3 and can be as high as 4.6% in people older than 75 years 4 of age [27]. Severe AS can also cause irreversible damage to 5
blood cells due to the augmented stress red cells and platelets 6
experience [26]. 7
Non-invasive Doppler echocardiography is the recom- 8
mended clinical standard to assess AS severity [3]. Guidelines 9
suggest to consider the peak transvalvular velocity, the mean 10
transvalvular pressure gradient calculated using a simplified 11
Bernoulli equation and the continuity-equation-derived aortic 12
valve area as the three primary haemodynamic parameters 13
for the assessment of AS severity [3]. Such commonly used 14
parameters for pressure drop estimation are based on strong 15
haemodynamic approximations and do not take into account 16
the interaction between valve, left ventricle and peripheral 17 vasculature ([13], [14], [22]). In fact, these indicators give 18 only partial information about the detrimental effects on the 19 left ventricle linked to the pressure loss caused by the mal- 20
functioning aortic valve [22]. 21
It has been reported that the transvalvular velocity peak 22
due to the reduction in cross-sectional area is not sufficient 23
to accurately estimate the additional pressure loss as it is 24
dependent on the flow conditions and a the threshold value on 25
velocity to evaluate stenosis severity would have to be then 26
patient dependent [3]. 27
The mean transvalvular pressure drop is calculated by means 28
of a modified or simplified Bernoulli law using only the 29
velocity averaged over the systolic phase in the vicinity of 30
the valve. This does not take into account the energy loss 31
due to the turbulent stresses brought about downstream from 32 the valve nor the energy loss across the stenosis orifice [28]. 33 As a result, for severe stenosis, the actual pressure drop is 34
underestimated when the simplified Bernoulli relation is used. 35
The degree of underestimation found in the literature varies 36
between 5% to 50% depending on the orifice area as well as on 37
the flow unsteadiness ([12], [15], [28]). In addition, it has been 38
shown that for mild to moderate stenosis with gradual distal 39
widening, the pressure drop localised at the valve calculated 40
with the simplified Bernoulli approximation overestimates the 41
actual pressure loss measured further downstream due to the 42
recovered pressure related to the decrease in velocity [3], [22]. 43
Other more complete models based on Bernoulli’s law have 44
been suggested and take into account energy losses across 45
the stenotic nozzle and downstream of the stenotic valve 46
This document is the accepted manuscript version of the following article:
Corso, P., Giannakopoulos, G., Gulan, U., Frouzakis, C. E., & Holzner, M. (2020). A novel estimation approach of pressure gradient and haemodynamic stresses as indicators of pathological aortic flow using subvoxel modelling. IEEE Transactions on Biomedical Engineering. https://doi.org/10.1109/
TBME.2020.3018173
along with the unsteady effects of flow acceleration during
1
systole([13], [28]). Many of these models use as independent
2
variable the effective orifice area (EOA) whose evaluation
3
based for instance on the continuity equation is a source of
4
inaccuracies and uncertainties [4]. Hatoum et al. [17] have
5
recently proposed a model that instead of the estimate of
6
EOA requires empirical constants which, however, might not
7
be suitable for a wide range of cases. The energy loss index
8
is a complementary indicator that accounts for irreversible
9
energy loss occurring in the aorta. The index is derived from
10
Bernoulli’s law and the integrated linear momentum equation
11
projected on the flow axis obtained under various assumptions,
12
i.e. uniform velocity profile, negligible acceleration term and
13
the estimation of EOA [13].
14
Finally, it is important to note that in addition to the assump-
15
tions made in order to formulate the aforementioned indicators,
16
measurement errors introduced in Doppler echocardiography
17
can be significant, due to the misalignment of the ultrasound
18
beam with respect to the flow direction resulting in a too high
19
intercept angle [3].
20
Although pressure loss evaluation from Doppler echocar-
21
diography can be easily implemented on patients in vivoand
22
can provide useful order of magnitude estimates of pressure
23
loss within a short acquisition time, in recent years, there has
24
been an effort to indirectly estimate pressure gradient and
25
energy loss by characterising the three-dimensional velocity
26
field usingin vivo4D flow magnetic resonance imaging (MRI)
27
and in vitro particle tracking velocimetry (PTV) [9], [14]–
28
[16], [19]. This has overcome some of the limitations of
29
Doppler echocardiography estimates. It is still based though
30
on several assumptions and the influence of limited spatial
31
resolution has never been investigated systematically. In fact,
32
4D flow magnetic resonance imaging (MRI) has recently been
33
employed to characterise turbulence in aortic flow and evaluate
34
irreversible pressure drop [15], [16]. Since pressure cannot be
35
directly measured from MRI measurements, Ha et al. [15]
36
assumed that the rate of work done by the pressure drop in
37
the stenosed aorta is equivalent to the total energy dissipated in
38
the flow with the latter estimated by the integral over the whole
39
aorta volume of the production of turbulent kinetic energy at
40
voxel sizes of 1 to 2 mm. A dedicated MR sequence has then
41
been developed and this sequence has been validated against
42
results obtained from large-eddy simulations to evaluate the six
43
components of the Reynolds stress tensor [15]. The product
44
of this reconstructed tensor with the shear-rate tensor obtained
45
from the MR-based velocity field gives the production of
46
turbulent kinetic energy. Moreover, it has been shown in [14]
47
and [16] that the latter can also be evaluated using a shear-
48
scaling argument requiring the use of a model constant and the
49
calculation of the mean shear rate and the diagonal elements
50
of Reynolds stress tensor. The limitations of such an approach
51
based on macroscopic energy balance lie in the applicability of
52
the aforementioned main assumptions to complex and bended
53
aorta geometries, to pulsatile flow conditions and to moving
54
boundary walls.
55
The approaches summarised above are based on several
56
assumptions to model pressure loss and focus on bulk pressure
57
loss, i.e. the loss over an extended control volume. To date,
58
there is no method available that is capable of predicting the 1 pressure gradient evolution across a given volume of interest 2
downstream of a stenosis. Furthermore, due to limited spatial 3
resolutions in flow measurement techniques such as particle 4
tracking velocimetry or 4D MRI, the evaluation of small-scale 5
stresses owing to turbulence in the flow still requires adequate 6
and validated models. 7
In this study, a direct numerical simulation of stenotic flow 8
using a spectral element flow solver is used as a benchmark 9
to establish a novel approach based on a rigorously derived 10
differential equation to evaluate pressure gradient along the 11
centreline of stenosed aortas. DNS data are thus used to 12
optimise models for subvoxel-scale stress estimation. In fact, 13
we test the hypothesis that the limited resolution inherent inin 14 vitroandin vivomeasurements, e.g. for 3D-PTV or 4D Flow 15 MRI, can be represented by a spatial filter and consequently 16 the under-resolved scales modelled with an LES-type model. 17
A systematic validation of the models for three voxel sizes is 18
performed. This new approach has been established based on 19
the characterisation of the flow in a stenosed aorta geometry 20
fromin silicoandin vitro three-dimensional particle tracking 21
velocimetry (3D-PTV) data. Finally, we also verify through 22
pressure catheter measurements the applicability of this new 23
approach to pulsatile flow conditions in a compliant aorta 24
model with a different stenosis geometry. 25
II. METHODOLOGY 26
Three methods are used to investigate stenotic flows. The 27
first investigated case concerns the flow downstream of a 28
severe stenosis of a tricuspid aortic valve in a rigid aorta 29
geometry (see Fig. 1). The combination of direct numerical 30
simulation,in silicoandin vitro particle tracking velocimetry 31
is applied to that flow case and allows us to develop and 32
optimise two models for subvoxel stress estimation as part of 33 our approach for pressure gradient evaluation. The validation is 34 performed for three different voxel sizes representative of the 35
resolution reached inin vivo4D Flow MR andin vitro3D-PTV 36
data, namely, 0.75, 1 and 1.5 mm. The second investigated case 37
deals with the flow after a concentric severe aortic stenosis in 38
the root of a compliant model of the same aorta geometry as 39
that of the first case. A pulsating flow rate is imposed upstream 40
from the stenosis and the approach based upon subvoxel-scale 41
models previously elaborated is validated using in vitro 3D- 42
PTV and pressure catheter measurements. 43
A. Direct Numerical Simulation (DNS) 44
The aorta geometry was obtained from an MRI scan and 45
has a diameter at the sinotubuluar junction of 25 mm. A 46
geometrical model of a stenosis was added to the sinotubular 47
extremity. This model represents a severe calcified stenosis of 48
a tricuspid aortic valve with a reduction in cross-sectional area 49
of about 80% ( Fig. 1). To generate the computational mesh, 50
the volume of the stenosed aorta lumen which is described 51
by a set of points on the external surface was reconstructed 52
by taking cross-sections perpendicular to the centreline and 53
using cubic splines to define the contours of each cross- 54
section of the lumen. The full geometry representation was 55
then completed by a sweeping operation of the reconstructed
1
cross-sections along defined path curves. The conforming
2
hexahedral curve-sided element mesh skeleton was generated
3
using the Trelis mesh generation toolkit (CSIMSOFT, 2018.
4
http://www.csimsoft.com/trelis.jsp). The domain was divided
5
into several blocks in order to allow for different mesh
6
refinements depending on the distance from the stenosis. The
7
final grid used consists of 92,208 conforming curve-sided
8
spectral elements within which the solution is expressed in
9
terms of order N=7 polynomials resulting in a total of about
10
47 million discretisation points.
11
Fig. 1. Geometry of the ascending aorta divided into different blocks to define different mesh cell sizes. The orifice geometry represents a severly stenosed tricuspid valve.
The evolution of the transitional flow is governed by the
12
incompressible continuity and momentum equations:
13
∇ •v= 0, (1)
∂v
∂t +v• ∇v=−∇p+ 1
Re∇ • ∇v, (2) with t, v, p, the nondimensional time, velocity and pressure
14
and the Reynolds number Re = U Dν = 2100 with U =
15
0.407m/s, the mean inflow velocity,D = 25mm, the inlet
16
diameter and ν = 4.85 10−6 m2/s, the kinematic viscosity.
17
The flow rate imposed duringin vitroPTV experiments is then
18
of 12 L/min.This flow rate corresponds to a typical systolic
19
flow rate pumped by an adult heart. The Reynolds number
20
calculated at the orifice Reo = vjetνdjet = 3800 with the
21
mean jet velocityvjet= 1.65m/sand the jet diameter at the
22
orifice djet = 11mm. The reference length and velocity for
23
nondimensionalisation are the inlet radius and twice the mean
24
inflow velocity, respectively. The simulations were carried out
25
using the open source incompressible flow spectral element
26
solver NEK5000 [21]. The solution is expressed in terms of
27
N-th order tensor-product Lagrange polynomials based on the
28
Gauss-Lobatto-Legendre quadrature points. The third-order
29
temporal integration is based on a semi-implicit formulation
30
using Backward Differentiation Formulas (BDF), where the
31
convective terms are treated explicitly and the viscous term
32
is treated implicitly [11]. A time-varying Dirichlet boundary
33
condition (BC) is imposed at the inflow, while the outflow
34
BC is similar to the one found in [29]. In order to ensure that 1 the flow characteristics at the outflow boundary are always 2
pointing outwards, the flow is forced to leave the domain by 3
imparting a positive velocity divergence in the last elements 4
attached to the boundary. Since the geometry of the domain 5
used for DNS is the exact same geometry as the one provided 6
for the manufacturing of the stenosed aorta model made of 7
rigid silicone and used in the experiments with the exception 8
that the domain is cut before the aortic arch and as the flow is 9
not pulsatile and no backflow at the level of the outflow bound- 10
ary is noticed during the experiments, this imposed outflow 11
BC is realistic and in line with the experimental conditions. 12
The non-trivial time-varying inflow velocity is obtained from 13
the measured profiles of each velocity component from in 14 vitro 3D-PTV experiments. Although the flow rate imposed 15 by the pump during the measurements is constant, it was 16 observed that the profiles were not perfectly steady in time. 17
In addition, the shape of these profiles is not axisymmetric 18
unlike the ones predicted by Poiseuille or Womersley solutions 19
and commonly used as inflow velocity profiles for simulations 20
[20], [29].We have noted that this unsteady behaviour caused 21
important differences in the flow development downstream of 22
the stenosis. In particular, the instability of the jet issuing 23
from the orifice varied significantly between a steady inflow 24
boundary condition and the time varying one (see video XY in 25
supplementary materials).For this reason, a proper orthogonal 26
decomposition (POD) using the method of snapshots for each 27
velocity component in the inflow plane was performed [5]. The 28
first 100 snapshots of this decomposition were used for the 29 reconstruction of the profiles. The number of retained modes 30 was chosen such that the sum of the 100 eigenvalues represents 31 95% of the total value over the full set of eigenvalues. The 32
reconstructed velocity profile is then expressed as: 33
viin(x, y, t) =
100
X
j=1
aj(t)Mi,j(x, y), (3)
where 34
Mi,j(x, y) =vimeas Vj
1
pλj, (4) is the jth mode for the ith velocity component and 35
aj(t) =p
λjVj, (5)
is the time-dependent coefficient multiplying each mode, while 36
λj and Vj are respectively the jth eigenvalue and the jth 37
eigenvector of the autocorrelation matrixRdefined as: 38
R=vTimeasvimeas ∈RNsnap×RNsnap. (6) The modes extracted from the POD are fitted to an analytical 39
function so that they can be imposed as boundary condition in 40
NEK5000. Since the inflow cross-section is circular, a radial
1
basis function (RBF) was chosen with the following form:
2
Mˆ =pkX
l
X
γ
exp(−r2+ 2rl(sinθsinγ+ cosθcosγ)−2l2) + exp(−r2−2rl(sinθsinγ+ cosθcosγ)−2l2) + exp(−r2+ 2rl(sinθsinγ−cosθcosγ)−2l2) + exp(−r2+ 2rl(−sinθsinγ+ cosθcosγ)−2l2) + exp(−r2),
l∈ {0.25,0.5,0.75}andγ∈ {0,π 6,π
4,π 3},
(7) with r the radius and θ the azimuthal angle of the inflow
3
cross-section. Through a fitting procedure, the coefficients
4
pk are determined such that the squared difference between
5
the value given by the RBF and the expected value of the
6
considered mode is minimised.The streamwise velocity profile
7
is constrained such that the integral of this component over
8
the inflow cross-section is equal to the flow rate measured
9
during the experiments. The flux for the other two velocity
10
components is constrained to be null. The shape of the RBF
11
that fits the reconstructed modes is shown at four time instants
12
in Fig. 2 (c), (d), (e), (f).
13
The smallest resolved scale in our DNS is about 25 microns.
14
The minimum Kolmogorov length scaleηalong the centreline
15
of the ascending aorta is equal to 75 microns and is located
16
at a distance from the stenotic orifice of about half the radius
17
at the sinotubular junction.
18
Within the scope of this study, we mimic the effect of
19
limited spatial resolution for in vitro and in vivo flow mea-
20
surements by spatially filtering the DNS data. Thereby, the
21
DNS velocity fields is filtered using an isotropic Gaussian
22
filter as the one commonly used for large-eddy simulation (see
23
Equations 11, 10, 21) [25]. The cut-off length∆ was chosen
24
to match the resolution usually achieved in in vitro PTV or
25
4D fow MRI measurements, namely 1.5 mm, 1 mm and 750
26
µm (see Fig. 3 (b)).
27
B. In silico 3D-PTV
28
The velocity along trajectory points obtained from three-
29
dimensional particle tracking velocimetry with perfectly pas-
30
sive tracer particles is extracted from the DNS velocity field.
31
The seeding density of tracer particles is chosen such that three
32
mean inter-particle distancesδmatch the previously mentioned
33
cut-off lengths ∆ (cf. Fig. 3 (c)). Thus, a seeding density
34
of about 250, 1000 and 2300 particles per cubic centimetre
35
results in a mean inter-particle distance of respectively 1.5, 1
36
and 0.75 millimetre (see video in supplementary material). In
37
order to calculate the relevant quantities for analysis from the
38
trajectory-based velocity, an interpolation of the velocity field
39
is performed using a Gaussian convolution with a window
40
width equal to 1.7δ [1]. In this work, the centreline of the
41
aorta geometry is the location where the three-dimensional
42
flow information is taken (cf. Fig. 3 (e)). In fact, pressure
43
catheter and Doppler measurements are usually made along
44
this centreline for the assessment of stenosis severity or for
45
the characterisation of the artificial valve performance. The
46
centreline can be determined as the weighted shortest path 1 between the two extremities of the ascending part of the aorta. 2
It is computed by minimising the line integral of the inverse 3
of the maximal inscribed sphere radius calculated at each 4
vertex of the Vorono¨ı diagram generated from the points on 5
the surface boundary of the ascending aorta [24]. 6
C. In vitro 3D-PTV 7
In vitrothree-dimensional particle tracking velocimetry ex- 8
periments are carried out in two anatomically accurate silicone 9
replicas of two stenosed aorta geometries. Details on the 10
apparatus and material used to make the measurements can be 11 found in [9] and [14]. The fluid used in both cases is a mixture 12 of water (48% by weight), glycerine (37% by weight) and 13 sodium chloride (15% by weight) that behaves as a Newtonian 14
fluid with kinematic viscosityν= 4.85 10−6m2/sand density 15
ρ= 1200kg/m3 and a refractive index of 1.41. Images are 16
recorded with a CMOS camera at a frame rate of 2 kHz and the 17
time interval between two instantaneous velocity fields after 18
the interpolation based on Gaussian convolution is around 3 19
ms. The geometry of the first rigid silicone model is the one 20
shown in Fig. 3 (f) and the one used in the direct numerical 21
simulation. The same nondimensionalisation of the fields of 22
interest is performed and a Gaussian convolution similar to that 23
applied to thein silicotrajectory point velocity is carried out. 24
As mentioned above, the velocities measured upstream of the 25
stenosis are used to impose the inflow boundary condition for 26 DNS. To test the applicability to a wider range of cases of the 27 new models proposed for subvoxel stress and pressure gradient 28 estimation,in vitromeasurements under pulsatile conditions in 29
the second aorta replica representing the same aorta geometry 30
made of flexible silicone are performed. The stenosis geometry 31
for this case is concentric and the reduction in cross-section 32
area is equivalent to that of the geometry presented in Fig. 1. 33
In vitro pressure measurements are also performed at various 34
locations using a sensor on a catheter. For further details on 35
the experimental setup and conditions, the reader is referred 36
to [14]. The interpolation method used to express the velocity 37
onto a Cartesian grid is similar to that of the previousin silico 38
andin vitroPTV case except that the time-averaging operation 39
is replaced by a phase-averaging one. 40
III. THEORY 41
The theoretical part of this paper starts with a rigorous 42 derivation of a differential equation for pressure gradient eval- 43 uation based on the Navier-Stokes equation and a description 44
of the filtering and ensemble-averaging operations applied to 45
it. Then, closure models for the subvoxel-scale stress term are 46
introduced. The final one-dimensional form of the pressure 47
equation is presented and the hypotheses made to formulate 48
the models are described. 49
A. Differential Equation for Pressure Gradient Evaluation 50 Since flow measurements resolved down to the Kolmogorov 51
length scale remain unavailable and since in clinical practice, 52
ensemble-averaged quantities are considered, we will assume 53
Fig. 2. (a) Ascending aorta geometry discretised with 47.2 106points. (b) Spectral element skeleton consisting of 92,208 conforming hexahedral elements.
Fitted inflow velocity profile retrieved from proper orthogonal decomposition (POD) of the velocity data obtained from in vitro PTV measurements at (c) t = 0 s, (d) t = 0.15 s, (e) t = 0.3, (f) t = 0.6 s.
that the behaviour of the flow quantities obtained with these
1
under-resolved flow measuring techniques is described by the
2
ensemble-averaged and spatially filtered Navier-Stokes vector
3
equation expressed as follows:
4
∇p˜
|{z}
Pressure gradient
= −˜v• ∇˜v
| {z }
Advective term
+ 1
Re∇ • ∇˜v
| {z }
Viscous term
−∇ •τ˜r
| {z }
Subvoxel-scale stress divergence
−∂v˜
∂t
| {z }
Acceleration term
(8)
where tilde refers to a filtered quantity:
5
φ(x, t) =˜ Z
Ω
φ(y, t)G(|x−y|)dy
= X
k∈Z3
φ(k, t) ˆˆ G(|k|) exp(ikx). (9)
The filter kernel G is Gaussian and defined in physical and
6
spectral spaces in a similar way to the one used for large-
7
eddy simulations [25]:
8
G(|x−y|) =
√12
√2π∆exp
−6|x−y|2
∆2
, (10)
G(|k|) = expˆ
−|k|2∆2 24
, (11)
with∆, the cut-off length or voxel size.
9
Ensemble-averaging is performed in the following way:
10
φ= 1 Nt
Nt
X
k=1
φinst(x, t+k∆t) (12) withNt, the number of time instances and∆tthe time interval 1
between those instances. In the case of constant flow rate, 2
Nt is equal to 450 and ∆t = 3 ms. Under pulsating inflow 3
conditions,Ntis the number of cycles considered which is 60 4
and∆t = 1 s is the period of the pulse. 5 By virtue of Reynolds decomposition, the advective term of 6
Eq. 8 is decomposed and the Reynolds stress term appears. The 7
ensemble-averaging operation is a Reynolds operator which is 8
commutative with time and space-derivative. 9
∇˜p=
−˜v• ∇˜v+ 1
Re∇ • ∇˜v− ∇ •R˜
−
∇ •˜τr
−∂v˜
∂t,
(13)
withR˜ the filtered Reynolds stress tensor: 10
R˜ij= ˜v0iv˜0j= ˜viv˜j−v˜iv˜j, (14) andv˜0 the filtered fluctuating velocity. 11
The penultimate term on the right-hand side of Eq. 13 repre- 12
sents the divergence of the so-called tensor of subvoxel-scale 13
(SVS) stress. This SVS stress is expressed in the following 14
way: 15
˜
τijr =^vijvij−˜vij˜vij. (15)
B. Model Derivation 16
Since the direct and rigorous computation of SVS stress is 17
only possible by filtering fully resolved DNS velocity fields, a 18
closure model based on resolved flow quantities is needed to 19
Fig. 3. Illustration of the three methods of investigation (a) DNS of the flow downstream of a severe stenosis geometry (method 1). Linear integral convolution visualisation of instantaneous velocity field coloured with the magnitude of pressure gradient in three planes shows flow patterns complemented with iso- surfaces of positive Q values highlighting coherent vortical structures. (b) Time-averaged flow features from DNS fields filter at three cut-off lengths (∆ = 1.5mm, 1 mm and 750µm) brought out by planes with linear integral convolution of velocity field and iso-surfaces of positive Q values. (c) Passive tracers advected using DNS velocity field (method 2). Three seeding densities were used such that the mean inter-particle distanceδis equal to the filter cut-off lengths of the filtered DNS. (d) Flow features as in (b) obtained from interpolation of velocity on the simulated trajectory points. The size of the voxel used for the interpolation is equal toδ. (e) Streamlines generated from the interpolated and time-averagedin silicovelocity field and reconstructed centreline of the ascending part of the considered aorta geometry. (f) Illuminated tracer particles duringin vitroexperiments in the rigid silicone model of the stenosed aorta (method 3). (g) Streamlines generated from the interpolatedin vitrovelocity field in the first half of the ascending aorta and centreline.
evaluate this quantity from e.g. in silico or in vitro 3D-PTV
1
fields.
2
The first model is based on the Smagorinsky-Lilly model
3
and uses the formulation suggested in [9]. In non-dimensional
4
form, it can be expressed as:
5
||˜τr||=Cτ,Smag r 1
Re
∆˜ q
||S˜inst||3, (16)
∇ •τ˜r=C∇•τ,Smag r 1
Re
∆˜
∇ q
||S˜inst||3
, (17) with ∆ = ∆/R, the non-dimensional voxel size and˜ S˜inst,
6
the filtered rate-of-strain tensor:
7
S˜inst=1
2 ∇˜v+∇˜vT
. (18)
A new model inspired by the mathematical development
8
presented in Peng and Davidson [23] is also proposed in this
9
work. A reconstruction series of the product of velocity vectors
10
filtered with an isotropic Gaussian kernel is first established.
11
This series is also referred to as the Leonard expansion [23].
12
The first term of this series is the matrix product of the velocity
13
gradient tensor with its transpose and is decomposed using the
14
resolved rate-of-strain and rate-of-rotation tensors (see Equa-
15
tions 18, 21, 22). It is also considered in this study to estimate
16
the fourth-order second term of the reconstruction series based
17
on the viscous dissipation rate tensor, whereas higher-order
18
terms are neglected. The tensor of viscous dissipation rate of
19
energy is evaluated by assuming that it is proportional to the
20
Hadamard product of the resolved time-averaged rate-of-strain
21
tensor and the Reynolds stress tensor, similar to the assumption
22
made in Ha et al.[15]. Therefore, the Frobenius norm of the
23
subvoxel-scale stress tensor and the last term on the right-hand
24
side of Eq. 13 are given by:
25
||˜τr||=
Cτ,1∆˜2
||˜Sinst˜Sinst||+Cτ,2∆˜2||Ω˜instΩ˜inst||
+Cτ,3∆˜2||S˜inst˜SinstΩ˜instΩ˜inst
||˜Sinst||2 ||
+Cτ,4Re∆˜2||˜SinstR˜
| {z }
P˜
||,
(19)
∇ •˜τr=
Cl,1∆˜2
∇ •˜SinstS˜inst+Cl,2∆˜2∇ •Ω˜instΩ˜inst
+Cl,3∆˜2∇ • ˜Sinst˜SinstΩ˜instΩ˜inst
||S˜inst||2
!
+CεRe∆˜2∇ •
˜SinstR˜ ,
(20) withΩ˜inst, the resolved rotation-rate tensor:
26
Ω˜inst= 1
2 ∇˜v− ∇˜vT
, (21)
and Equations 19 and 20 stem from the decomposition of 1 the product of the velocity gradient tensor with its transpose: 2
∇˜v∇˜vT =S˜inst˜Sinst+Ω˜instΩ˜inst +˜Sinst˜SinstΩ˜instΩ˜inst
||S˜inst||2 .
(22)
In addition, in line with the objective of this study, the 3
pressure gradient vector needs to be projected along the 4
aorta centreline. The differential equation for pressure gradient 5
evaluation can be written in its final one-dimensional form: 6
∇p˜•t=
−˜v• ∇˜v+ 1
Re∇ • ∇˜v
− ∇ •R˜
•t
−
∇ •˜τr
•t−∂v˜
∂t •t.
(23)
witht the tangent vector to the centreline. The last term on 7
the right-hand side of Eq. 23 is equal to zero in the case of 8
constant imposed flow rate. 9
The approach followed in this study consists of optimising 10
the parameters of Equations 23, 16, 17, 19, 20 for different 11 voxel sizes by calculating the various termsin the right-hand 12 side of these equations from in silico PTV velocity fields 13 and projecting them onto the tangent vector of the centreline. 14
Pressure gradient fields from filtered DNS as well as subfilter- 15
scale stress retrieved from DNS data are used as reference 16
values. Thus, the left-hand side terms of Equations 16, 19, 23 17
are obtained from filtered DNS fields whereas the terms on 18
the right-hand side of these equations are all evaluated from 19
in silico or in vitro PTV. The parameters Ci appearing in 20
Equations 16, 17 19, 20 are obtained by a linear least-square 21
optimisation using an interior-point method [2] and performed 22
over a set of around 4500 points located in a tube centred 23
on the centreline and under constraints on the integral of the 24
various quantities along the centreline. This problem can be 25
formulated in the following way: 26
min
Ci f(Ci) = min
Ci ||M Ci−d||22
= min
Ci
1
2CiTHCi+rTCi such thatAeqCi=beq.
(24)
with f(Ci), the objective function; M, the matrix with the 27
terms on the right-hand side of Equations 16, 17, 19, 20, 23; 28 d, the vector with the reference value of pressure gradient 29
and subvoxel-scale stress obtained from DNS; H, the Hessian 30 matrix2MTM;r=−2MTd;Aeq, the matrix with the values 31
of the integral along the centreline of the different terms; 32
beq, the reference pressure difference and the integral over the 33
centreline of the reference SVS stress. To measure how close 34
the computed minimum is to the optimal value as part of the 35
interior-point algorithm, the first-order optimality is tracked. A 36
relative tolerance of10−14for this measure is set as a stopping 37
criterion for the iterations. The first-order optimality measure 38
when constraints are present is based on the Karush-Kuhn- 39
Tucker conditions, which can be expressed as: 40
||∇CiL(Ci, λ)||∞=||∇f(Ci) +ATeqλ||∞ (25)
with λ, the Lagrange multiplier vector from the equality
1
constraint on the integral values.
2
In order to compare the performance of the two proposed
3
models for the three considered resolutions, the coefficient
4
of determination is calculated. As defined in Eq 26, this
5
coefficient represents the percentage of the variance of the
6
observed data (i.e. the pressure gradient and the SVS stress)
7
explained by the model. The computation of such a coefficient
8
is performed based on a set of about 4500 points distributed in
9
a tube centred on the centreline of the aorta and with a radius
10
of approximately 0.8 times the aorta radius. The reason for
11
choosing this volume is to have a sufficient number of data
12
points for a robust optimisation of parameters general enough
13
to estimate SVS stress and pressure gradient in the bulk of
14
the flow. Indeed, a convergence study on the R2 values was
15
performed to make sure that the number of sampling points
16
for the optimisation is sufficient. We will show in the results
17
section below the estimated pressure gradient and SVS stress
18
along the centreline, which is a subset of the whole set of
19
points considered for optimisation.
20
R2= 1− P
i
φi−Bˆi
2
P
i φi−φi2, (26)
Bˆiis the expected regression line with slope equal to 1,φiare
21
the values of pressure gradient or SVS stress evaluated with
22
the models at the sampling points and φi is the of average
23
value of the reference pressure gradient or SVS stress.
24
Furthermore, in order to assess the sensitivity of the two
25
models at the three resolutions to a variation in the parameter
26
value or error in the evaluated terms, we study the variation of
27
the coefficient of determination as a function of an imposed
28
error ranging from -40% to + 40% which is a conservative
29
choice because it is larger than typical experimental errors in
30
velocity data and related terms of Eq. 23 ([9], [14], [19]).
31
Since the models for pressure gradient estimation and the new
32
model for SVS stress evaluation use more than one optimised
33
parameter, a combined R2 is calculated by computing the
34
convolution product of the 6th degree polynomials fitted to
35
the R2 distribution depending on the imposed error for each
36
term of the models (cf. Equations 23, 16, 17, 19, 20).
37
Finally, our approach is based upon the assumption that the
38
velocity field obtained from the Gaussian filtering operation
39
applied to the DNS velocity is consistent with the velocity
40
field interpolated on a structured Cartesian grid drawn from the
41
Gaussian convolution applied to the trajectory point velocities.
42
IV. RESULTS ANDDISCUSSION
43
This section presents the results of the model constant opti-
44
misation and compares the estimated pressure and SVS stress
45
with the reference DNS-based values to draw conclusions
46
on the validity of models in relation to the three considered
47
voxel sizes. It also provides a physical interpretation for the
48
new model and discusses the applicability of these models
49
to in vitro PTV measurements. It then compares pressure
50
drop estimates obtained with other pressure models from the
51
literature that are used in clinical practice. Finally, the validity
52
of the new approach is tested for more realistic conditions and 1 the relevance for clinical application of such a novel approach 2
is highlighted by comparing the discrepancies from the direct 3
pressure measurement between the proposed models and the 4
method proposed by Haet al.([15], [16]). 5
A. Model Optimisation, Performance and Physical Interpre- 6
tation 7
Tables I and II provide the values of the optimised coeffi- 8
cients. 9
The performance of the modelled projected pressure gra- 10 dient and subvoxel-scale (SVS) stresses as compared to the 11 reference values is assessed by considering the coefficient 12 of determination R2 given in Tables I and II calculated 13
based on a sample of 4500 points. The optimisation of the 14
above parameters is performed for the Equations 17, 20, 23 15
allowing for the evaluation of pressure gradient but also for 16
the Equations 16, 19 representing the norm of the SVS stress 17
tensor. The usefulness of the SVS stress models comes from 18
the blood damage models, which are employed to quantify the 19
detrimental effect blood flow exerts on the blood cells along 20
their trajectory [9]. 21
As far as the Smagorinsky-based models are concerned, we 22
note in the first row of Tables I and II that the prediction 23
quality of both the pressure gradient and the SVS stress 24
reconstructed fromin silico PTV is independent of the voxel 25
size. The coefficient of determinationR2representing the pro- 26
portion (as a percentage of the variance of the reference data) 27
of the in silico PTV estimate that agrees with the reference 28
data is around 85% for the SVS stress estimate and 60% for 29
the projected pressure gradient estimation. This trend suggests 30
that irrespective of the voxel size and for optimised constants 31 C∇•τ,Smag and Cτ,Smag that do not strongly depend on the 32 voxel size, these models would provide a good estimation 33
especially for the SVS stress. 34
It is worth noting that the lowest resolution considered 35
in this study corresponds to the one that can commonly be 36
achieved with current in vitroand in vivo flow measurement 37
techniques such as 3D-PTV or 4D-flow MRI [14], [15], [19]. 38 With regard to the new model for the evaluation of pressure 39
gradient from in silico PTV, the prediction quality is also 40
independent of the the decreasing voxel size as evidenced by 41
the values of R2 in Table I. If we compare the performance 42
of this new model with the Smagorinsky-based one, we notice 43
that the coefficient of determination is 25% larger. Overall, the 44 improvement in terms of fitting quality made with the new 45
model is substantial. 46
The prediction quality for in silico SVS stress using the 47
new model proposed in this study is reported in the first row 48
of Table II. The coefficient of determination is 14% smaller 49
at the two higher resolutions than the one at the voxel size of 50 0.12. This suggests that this new model performs better when 51
the resolution is coarser. 52
The physical interpretation of the new model for the eval- 53
uation of pressure gradient originates from the fact that in 54
turbulent shear flows, many of the important kinematic and 55
dynamical properties can be understood by describing the 56
TABLE I
PARAMETERS OPTIMISED FOR THE TWO MODELS CONSIDERED FOR THE EVALUATION OF PRESSURE GRADIENT PROJECTED ONTO THE
CENTRELINE.
∇p•t ∆˜ Cl,1 Cl,2 Cl,3 Cε C∇•τ,Smag R2 ||∇CiL(Ci, λ)|| ×10−14 λ Smagorinsky
0.12
n/a n/a n/a n/a
13.25 0.64
n/a n/a
0.08 12.03 0.61
0.06 8.74 0.65
New model
0.12 1.26 -0.63 0.028 0.013 n/a
0.77 0.96 -0.014
0.08 1.96 2.14 0.095 0.094 0.75 2.45 -0.305
0.06 2.07 0.5 -0.0064 0.196 0.78 0.116 0.0126
n/a = not applicable.
TABLE II
OPTIMISED PARAMETERS FOR THE MODELS ALLOWING TO CALCULATE THE MAGNITUDE OF SUBVOXEL-SCALE STRESS ALONG THE CENTRELINE
FROM3D-PTV.
||τr|| ∆˜ Cτ,1 Cτ,2 Cτ,3 Cτ,4 Cτ,Smag R2 ||∇CiL(Ci, λ)|| ×10−14 λ Smagorinsky
0.12
n/a n/a n/a n/a
21.71 0.87
n/a n/a
0.08 21.66 0.84
0.06 23.59 0.82
New model
0.12 0.00 6.31 -3.2 -0.018 n/a
0.89 0.22 -0.205
0.08 0.00 4.86 1.13 0.03 0.76 5.86 6.68
0.06 0.00 5.52 -1.52 0.1 0.77 2.84 6.62
flows in terms of individual events, namely the combination
1
of coherent eddies, convergence and streaming zones [18].
2
In that respect, qualitative and mathematical criteria based
3
upon threshold values for the second invariant of the velocity
4
deformation tensor and the pressure in the flow have been
5
proposed in Huntet al.[18]. This representation can be useful
6
since in this study, ensemble-averaged and under-resolved
7
features of a transitional-to-turbulent flow can be characterised.
8
The ensemble-averaged description of the flow at a fixed
9
voxel size is then described through a combination of events
10
representative of the fundamental nature of the underlying
11
turbulence. Following this train of thought, because the new
12
model for pressure gradient evaluation links pressure gradient
13
and velocity deformation tensor through an equation, we sense
14
that the optimisation of the different parameters could lead
15
to the definition of unequivocal criteria that would allow to
16
highlight the interplay between individual patterns in the type
17
of shear-induced turbulent flows encountered in the aorta.
18
Thereby, each term of Equations 13 and 20 can be interpreted
19
in the following manner:
20
• −˜v • ∇˜v represents the contribution to the pressure
21
gradient in the flow due to streaming zones. As stated
22
by Hunt et al.[18], in streaming zones, the flow is fast,
23
not very curved and not diverging or converging strongly.
24
• −∇ •R˜ represents the contribution to the pressure gra-
25
dient due to momentum transport from Reynolds shear
26
stresses.
27
• −∆˜2∇ • ˜Sinst˜Sinst represents the contribution to the
28
pressure gradient due to irrotational straining motion
29
in the flow also known as convergence zones (TCl,1).
30
These zones are characterised by strong convergence and
31
divergence of streamlines.
32
• −∆˜2∇ •Ω˜instΩ˜inst represents the contribution to the
33
pressure gradient due to the action of coherent eddies
34
in the flow (TCl,2) defined as strong swirling zones with
35
vorticity.
36
• −Re∆˜2∇ •
˜SinstR˜
represents the contribution to
37
the pressure gradient due to turbulent energy dissipation
38
(TCε).
39
TABLE III
CONTRIBUTION OF EACH TERM OF THE EQUATION FOR PRESSURE GRADIENT ESTIMATION FOR THE TWO VARIANTS OF THE NEW MODEL. THE CONTRIBUTION IS CALCULATED FOR EACH VOXEL SIZE FROM THE RATIO BETWEEN THE INTEGRAL OF EACH TERM OVER THE CENTRELINE
AND THE VALUE OF RECOVERED PRESSURE.
R lTi•tdl
∆precovered
∆˜ Tv•∇˜˜ v T∇•R˜ TCl,1 TCl,2 TCl,3 TCε
New model
0.12 -2.46 1.61 1.93 0.35 -0.01 -0.43 0.08 -2.58 1.75 3.58 -0.39 -0.03 -1.36
0.06 -2.6 1.92 3.77 -0.1 -0.002 -2.09
Therefore, the new model weighs each of the last three terms 1
above such that the combined action of the physical flow 2
characteristics best matches the reference pressure gradient 3
value. The contribution of each weighted term is presented in 4
Table III. For this, each term is integrated over the centreline 5
and divided by the value of recovered pressure in the ascending 6
part of the considered stenosed aorta. 7
The term representing the influence of the turbulent straining 8
motion zones on the pressure gradient prevails for the three 9
voxel sizes and tends to counteract the action of the streaming 10
advection (negative values) especially close to the stenosis. 11
The absolute value of the relative contribution due to turbulent 12
energy dissipation term (TCε) increases with the diminishing 13
voxel size whereas the streaming zone term and the Reynolds 14
stress term remain rather constant for the three voxel sizes. 15
B. Sensitivity Analysis of The Models 16
Figure. 4 presents the variation in the coefficient of deter- 17 mination for the two models and for the estimation of both the 18 pressure gradient and SVS stress at the three resolutions. This 19
analysis is useful to understand the impact of measurement 20
errors or change in the parameter value on the predicted quan- 21
tities of interest. The noticeable trend is the more important 22
sensitivity of both models at the coarsest tested resolution 23
whether for the pressure gradient or SVS stress estimates. In 24
particular, as seen in Fig. 4 (a), the Smagorinsky-based model 25
for the estimation of SVS stress magnitude is very sensitive to 26
an over-estimation of its only term which is dependent on the 27
calculation of the strain rate tensor magnitude. That being said, 28
the new model for SVS stress estimation is much less sensitive 29
even at the coarsest voxel size and for an error reaching 40% 30
of the optimal state, the decrease in theR2is limited to a value 31 of about 50%. It is worth noting that unlike the Smagorinsky 32 model at the larger voxel size, the new model for the evaluation 33
of SVS stress at the intermediate and smaller voxel size is 34
more robust to a general over-estimation than under-estimation 35
of the various terms. With regard to the sensitivity of the 36
two models for the pressure gradient estimation, the quality 37
of prediction with the new model at the smaller voxel size 38
is limited to a decrease in R2 value of 20% when a 40% 39
error is introduced. The Smagorinsky-based is more sensitive 40
to errors or variation in its parameter value than the new 41
model except at the lower resolution where the new model 42
reaches negative R2 values when a 30% over-estimation or 43
under-estimation of its composing terms is introduced. Finally, 44
in order to identify the term with the strongest influence on 45
