On the numerical treatment of viscous and convective effects in relative pressure reconstruction methods.

The mechanism of many cardiovascular diseases can be understood by studying the pressure distribution in blood vessels. Direct pressure measurements, however, require invasive probing and provide only single-point data. Alternatively, relative pressure fields can be reconstructed from imaging-based velocity measurements by considering viscous and inertial forces. Both contributions can be potential troublemakers in pressure reconstruction: the former due to its higher-order derivatives, and the latter because of the quadratic nonlinearity in the convective acceleration. Viscous and convective terms can be treated in various forms, which, although equivalent for ideal measurements, can perform differently in practice. In fact, multiple versions are often used in literature, with no apparent consensus on the more suitable variants. In this context, the present work investigates the performance of different versions of relative pressure estimators. For viscous effects, in particular, two new modified estimators are presented to circumvent second-order differentiation without requiring surface integrals. In-silico and in-vitro data in the typical regime of cerebrovascular flows are considered, allowing a systematic noise sensitivity study. Convective terms are shown to be the main source of error, even for flows with pronounced viscous component. Moreover, the conservation (often integrated) form of convection exhibits higher noise sensitivity than the standard convective description, in all three families of estimators considered here. For the classical pressure Poisson estimator, the present modified version of the viscous term tends to yield better accuracy than the (recently introduced) integrated form, although this effect is in most cases negligible when compared to convection-related errors.

INTRODUCTION

Blood pressure is not only an important clinical marker, but also a determining factor in the onset and evolution of cardiovascular diseases. This is true not only for the absolute magnitudes of the pressure—for example, measured through catheterization—but also for how the pressure is distributed within a vessel relatively to a reference point in its lumen. Relative pressure can be understood as a flow's response to inertial effects, viscous resistance and volume forces , ; yet, differently from flow velocities, it cannot be visualized or measured via imaging techniques. Thus, estimating relative arterial pressure from blood‐flow velocities is a relevant task in inverse hemodynamics. As matter of fact, recent studies comparing computed and measured relative pressure show that reconstruction methods can be reliable alternatives to invasive catheter‐based techniques. , , , , , Furthermore, pressure reconstruction is not only relevant in biomedical, but also in general engineering applications using, for example, particle‐based velocimetry. , ,

Depending on the application, there are essentially two classes of relative pressure estimators: those capable of computing the full pressure field in space and time, and those providing mean pressure differences between selected planes throughout time. Pressure drop markers have a well‐established clinical importance as a metric for stenotic severity, and can also indicate the risk of growth and rupture in aneurysms. , , Full‐field relative pressure estimators, on the other hand, provide more detailed information with potential diagnostic/prognostic implications: the spatial distribution of intra‐arterial pressure, for instance, is correlated with where an infundibulum is likely to rupture. , It is important to remark, however, that a homogeneous pressure rise may also increase the risk of aneurysm rupture but can unfortunately not be captured by relative pressure estimators. This limitation stems from the fact that incompressible flow equations contain only the pressure gradient, so that determining absolute values would require either a reference pressure measurement (at one point at least), or some information on, for example, vessel wall compliance. Therefore, the pressure estimation methods addressed herein should be understood as relative estimators, an overview of which is found in the recent article by Bertoglio et al.

The most common application of pressure estimators is in determining pressure drops across stenoses, especially in the aorta. , , , , , , Because of the large flow velocities therein, the pressure gradient tends to be driven predominantly by convective effects. On the other hand, advances in imaging techniques in the past decade have enabled accurate 4D flow measurements also in smaller vessels, , , , , , , where viscous effects tend to be more significant. These effects may be even more pronounced in diseased arteries: for instance, typical flow structures in intracranial aneurysms resemble that of lid‐driven cavity flows with strong shear layers. In fact, relative pressure reconstruction in cerebral vessels has gained increased attention in the last few years, , , , although further clinical analysis on its impact is still required. Diffusive effects can play a major role not only in low Reynolds regimes, but also in turbulent flows where an eddy viscosity is introduced to model turbulent dissipation, and pressure estimators incorporating turbulent effects have been recently presented. , ,

Despite the mathematical simplicity of the viscous term in the (Newtonian) Navier–Stokes system, pressure estimation methods often neglect this term because of its higher‐order velocity derivatives. In the popular pressure Poisson estimator (PPE), for example, second‐order spatial derivatives are present even when using weak formulations. Therefore, retaining viscous effects normally requires higher‐order interpolation of the velocity data. Besides, second‐order differentiation can lead to large errors for noisy velocity measurements and should thus be avoided when possible. Pacheco and Steinbach modified the PPE to include the viscous contribution using only first‐order derivatives, but the required boundary integral could reduce accuracy in the presence of near‐wall data distortions and imaging artifacts. The Stokes estimator (STE) recently introduced by Švihlová et al. and revisited by Bertoglio et al. can account for viscous effects without second‐order derivatives or boundary operators, but at the cost of quadrupling the number of unknowns in comparison to the PPE. In spite of its higher cost, the Stokes estimator is among the simplest and most accurate techniques for reconstructing relative pressure fields from velocities.

While viscous forces can be troublemakers due to higher‐order differentiation, convective forces can also lead to increased noise sensitivity due to data squaring. The standard form of the convective term contains the product between the velocity field and its spatial gradient, which is why Bertoglio et al. proposed a modification of the Stokes estimator integrating this term by parts to get rid of velocity derivatives. Also proposed therein is the integral momentum relative pressure (IMRP) estimator where integration by parts is applied to reduce the order of differentiation of both viscous and convective terms. Yet, shortly after that, Marlevi et al. reported that the surface operator produced by integrating the convective term by parts can deteriorate the method's accuracy, and thus proposed the virtual work‐energy relative pressure (WERP) estimator keeping the convective term in its non‐integrated form. To the best of the author's knowledge, however, no actual quantitative comparison between IMRP and WERP estimators is yet available. Moreover, applications of both integrated and non‐integrated versions of the Stokes estimator can be found in the recent literature, and only a limited comparison has been reported to date.

Given the importance of both viscous and convective terms for the accuracy of relative pressure estimators, as well as the shortage of comparative studies assessing different ways of treating such terms, this work presents a systematic investigation considering three popular families of estimators in different variants. Classical and state‐of‐the‐art techniques are considered, as well as new modified treatments of the viscous term. For the PPE, a numerical technique is introduced to reconstruct viscous forces without resorting to boundary integrals or second‐order derivatives; for the WERP estimator, a simple modification based on the rotational form of the momentum equation is proposed to eliminate the viscous boundary term. Different test cases are carried out with increasing complexity, starting from two‐dimensional flow problems with synthetic measurements and then moving to more practical tests considering 4D magnetic resonance imaging (MRI) data. It is important to remark that when using other imaging modalities, such as ultrasound, a three‐dimensional velocity field might not be immediately available but can be reconstructed using appropriate techniques. , In that case, the pressure can be estimated either after reconstructing the velocity field, using for example the methods discussed herein, or reconstructed with the velocity. The reader is referred to the recent work of Galarce et al. for details on state‐of‐the‐art reconstruction techniques from Doppler ultrasound imaging.

PRELIMINARIES

An assumption common to most pressure estimators is that the flow velocity and the pressure satisfy the Navier–Stokes momentum equation: where is the fluid's density and is the viscous stress tensor, which under the Newtonian assumption is given by with denoting the dynamic viscosity. The pressure gradient is thus driven by three contributions: viscous, transient and convective forces. Both the viscous and the convective terms have the potential to increase noise sensitivity in pressure reconstruction: the former due to its second‐order derivatives, and the latter because of the “squaring” of the velocity data. The focus of the present work is thus placed on different treatments of these terms.

The conservation of mass for an incompressible flow reads simply which allows different simplifications of the stress divergence , as well as consistent modifications of the convective term. For the viscous term, the most popular version is the Laplacian form, with denoting the (vector) Laplace operator. Although measured hemodynamic velocities are in practice not divergence‐free—because of measurement and interpolation errors—, the simplification above is still consistent. Also noting that where denotes the curl operator, leads to the rotational variant: Although not so popular in flow simulations due to its somewhat unusual boundary conditions, this formulation is especially interesting for pressure reconstruction, as will be shown later on.

As for convection, the standard (often called “convective”) form can be recast into the so‐called “divergence form” or “conservation form”: since This variant arises, for example, when applying integration by parts to Stokes or virtual‐work‐based , estimators. Although all of these forms are equivalent when considering perfect (divergence‐free) measurements, they can perform very differently in the presence of noisy data.

The focus herein are techniques relying on Lagrangian interpolation of the velocity data, in particular using continuous, piecewise (tri‐)linear finite element spaces . The interpolated velocity field at time will be denoted by . A backward finite difference is used for the temporal discretization: where . Consider a spatial domain , or 3, with a boundary split into three non‐overlapping regions: , denoting inlet, outlet and vessel wall, respectively. In stenotic vessels, for instance, an important clinical measure is the mean pressure drop across the obstruction: When studying other conditions such as aneurysms and infundibula, there is also interest in estimating the actual relative pressure field. In that case, the vector‐valued momentum Equation (1) becomes an overdetermined system, since the velocity field is already known and one is left with equations to find a single scalar field. Different estimators arise depending on the method used for solving this inverse problem, as will be seen throughout this work.

CLASSICAL AND STATE‐OF‐THE ART ESTIMATORS

Pressure Poisson estimators (PPE)

The most classical technique able to fully quantify the pressure field is the so‐called pressure Poisson estimator (PPE). It results from the optimization problem to find the minimum of as shown by Song et al. Using a variational principle leads to the weak Poisson problem of finding such that with consistent Neumann boundary conditions already implied. Since this is a pure Neumann problem, the pressure at each time is defined up to an arbitrary additive constant. An additional constraint is thus required, for instance enforcing that A critical limitation of the standard PPE is the viscous term containing second‐order spatial derivatives of . Standard Lagrangian finite elements are only ‐continuous and therefore unable to account for higher‐order derivatives (regardless of the polynomial degree used in the interpolation!). Therefore, computing the viscous term requires either ‐continuous interpolation of the flow data or finite differences. To avoid large errors induced by second‐order differentiation, it is also common to simply neglect the viscous term, due to the typically small viscous contribution in stenotic flows.

To allow capturing viscous effects without computing second‐order derivatives and within a Lagrangian finite element framework, a least‐squares‐like reconstruction step is proposed here for the stress divergence. By projecting the viscous tensor onto a continuous finite element space before introducing it into the PPE, the resulting tensor will lead to a stress divergence which is piecewise constant, hence non‐vanishing. More precisely, first compute each component , with , such that with denoting the ‐th spatial component of , and considering the quantities evaluated at the desired time instant . Then, find such that Notice that, due to the symmetry of the stress tensor, only (rather than ) scalar projections are needed. In three dimensions, this corresponds to six scalar mass‐matrix systems—all with the same coefficient matrix—that must be solved at each time step. The added computational cost can be greatly mitigated by storing and re‐using the LDL decomposition of the mass matrix. This projected‐stress estimator (PPE) has the advantage of allowing non‐homogeneous viscosity, as arising from turbulent or non‐Newtonian hemodynamic models.

A technique sometimes employed to circumvent second‐order differentiation is integrating the viscous term by parts. This is however not immediately possible in pressure Poisson estimators, since the test functions already have first‐order derivatives in the weak form. Nonetheless, using the rotational description of the viscous forces (5) actually enables such a step. Notice that since in even for standard functions. This leads to the boundary vorticity or integrated pressure Poisson estimator (PPE) recently proposed by Pacheco and Steinbach, seeking such that where denotes the unit outward normal vector on . Although simple and computationally cheap, this estimator relies on a surface integral, a type of operator which is sometimes avoided in MRI‐based pressure reconstruction due to typically inaccurate near‐wall flow measurements. The performance of these two variants (PPE and PPE) will be compared in Section 4, including a third variant PPE using the divergence form (6) of the convective term.

Stokes estimators (STE)

The overdeterminancy of (1) can be handled by introducing a vector‐valued unknown and an additional scalar equation. In that context, and to circumvent the higher regularity requirements of the classical PPE, Švihlová et al. proposed perturbing the momentum equation with the Laplacian of an (artificial) incompressible velocity field . This leads to the Stokes problem of finding such that When using finite element methods to solve this problem, one can take advantage of the variational framework and use integration by parts to get rid of the second‐order derivatives in . In fact, it is also possible to eliminate the derivative in the convective term: the integrated Stokes estimator (STE) proposed by Bertoglio et al. seeks , with , such that for all , with , where is a generic stabilization term to allow equal‐order finite element spaces . Notice that the boundary terms coming from integration by parts vanish because . Of course, other forms of the viscous term can also be used, such as as proposed by Pacheco and Steinbach, or even the rotational form: in which case we will use the notation STE. Although there is no obvious a priori reason to pick one over the others – especially when has been treated with solenoidal filters/projections , , –, the “elasticity form” (19) has the advantage of allowing variable viscosity. Nevertheless, only the Laplacian (17) and rotational (20) forms will be considered herein.

Work‐energy balance estimators

When the goal is not to reconstruct the entire relative pressure field, but only to estimate an average pressure drop between two predefined planes, the effort in solving a boundary value problem for the pressure can actually be spared. Bertoglio et al. recently proposed using a virtual divergence‐free vector as weighting function for the momentum Equation (1): By taking such that and , integration by parts leads to: Furthermore, Therefore, defining the virtual flow rate and the ‐weighted pressure drop yields the integral momentum relative pressure (IMRP) estimator: with the virtual field typically pre‐computed as a Stokes flow. A variant of this method, the virtual work‐energy relative pressure (WERP) estimator does not rely on integration by parts of the convective term: In particular, setting recovers the original work‐energy (WERP) estimator by Donati et al. Their method is less accurate than its virtual relatives, not only because of the squaring of the velocity data, but also due to limiting assumptions (e.g., non‐bifurcating geometry, immovable vessel walls with purely tangential velocity); it can also be inaccurate for low flow rates and has a provenly higher bias than other estimators (refer to Bertoglio et al. for a detailed discussion).

As in the STE, in principle any equivalent form of the viscous term can be used here. However, for the virtual work‐energy estimators, the rotational form brings a potential advantage. In practical applications of the WERP estimator, the viscous boundary integral in Equation (27) (often denoted as ) is usually neglected, in order to mitigate errors related to inaccurate near‐wall measurements. , , , Although this is in practice a reasonable assumption, especially in typical cardiovascular flow regimes, it will be shown next how it is possible to construct a scheme where the boundary integral naturally vanishes. When applied to the virtual work‐energy equation, the rotational form of the stress divergence yields: Compare the boundary term in Equation (28) to the one in Equation (27). In the original WERP estimator, that term is proportional to the full trace , so that eliminating the boundary integral would require a condition—thereby also annihilating the pressure term on the left‐hand side of (27). Conversely, the rotational version (28) contains not the full but only the tangential trace , so that the boundary integral can be easily eliminated by selecting a virtual field with normal inflow and outflow conditions. For example, can be chosen as the solution of the Stokes problem to find satisfying with being a given scalar function on with zero value on and zero mean on (to guarantee global mass conservation for ). Then, it is simple to verify that the boundary integral in Equation (28) vanishes, since In that case, one gets the rotational (WERP) estimator: in which In particular, uniform inflow and outflow profiles can be used: This specific choice of leads to: that is, the weighted pressure drop is replaced by the standard mean pressure drop.

COMPUTATIONAL TESTS

In this section, various computational tests are carried out to evaluate the performance of the different estimators described so far. First, two‐dimensional problems with analytical solution are employed to assess the sensitivity of the different estimators with respect to noise. To isolate the influence of viscous and convective terms, the first and second examples are Navier–Stokes flows that also satisfy, respectively, the Stokes and Euler equations. Finally, 4D flow MRI data from an in‐vitro experiment will be used as a more practical test case. To assess noise sensitivity, the linearly interpolated velocity field is perturbed at each time step by an additive noise field with a normal distribution , so as to simulate real measurements. , The noise level is relative to the (spatiotemporal) peak velocity. Given the large number of possible versions for the estimators depending on how viscous and convective terms are treated, the notation used throughout this article is given in Table 1.

TABLE 1

Notation employed to describe the three classes of estimators (Poisson, Stokes and work‐energy) with respect to different treatments of viscous and convective effects

Abbreviations	Viscous terms	Convective terms
PPE	Neglected	Convective form
PPE_S	Stress divergence	Convective form
PPE_ω	Rotational form (integrated)	Convective form
PPE_div	Stress divergence	Divergence form
STE	Laplacian form (integrated)	Convective form
STE_int	Laplacian form (integrated)	Divergence form (integrated)
STE_int,ω	Rotational form (integrated)	Divergence form (integrated)
IMRP	Laplacian form (integrated)	Divergence form (integrated)
vWERP	Laplacian form (integrated)	Convective form
vWERP_ω	Rotational form (integrated)	Convective form

In the computations, first‐order Lagrangian finite elements are employed for the discretization of all quantities of interest. To solve the Stokes system required for the virtual work‐energy and Stokes estimators, a stabilization method is required to allow equal‐order discretization. For the two‐dimensional problems, which are solved in triangular meshes, the MINI element is used. , Since these elements are not defined for non‐simplicial meshes, in the 4D MRI case the stabilization method by Dohrmann and Bochev is employed — both techniques are absolutely stable, parameter‐free and simple to implement. All computations were done in MATLAB (MathWorks, Natick, MA) using its standard direct solvers.

Two‐dimensional Womersley flow

Given their known analytical solutions, two‐dimensional Womersley flows are a useful benchmark for pressure estimation methods. , The setup considers a pulsatile velocity given by in a straight channel , where and denotes the real part of a complex number. The exact pressure field is (up to a function of time only) so that the analytical expression for the pressure drop is The parameters used here are  m/s,  rad/s,  kg/m3,  mPa s, and  mm. Since the Womersley flow is not only a Navier–Stokes but also a Stokes flow, the idea with this test case is to at first drop the convective term in the estimators so as to compare only the different treatments of the viscous term. The triangular mesh has a uniform spatial resolution of , and the temporal resolution is set as , with s being the pulsation period. The estimators' sensitivity to noise is evaluated through the relative peak error: with 30 samples (cycles) computed per noise level. Two versions of each family of methods are considered: PPE and PPE, WERP and , and also the STE in standard (grad‐grad) and rotational (curl‐curl) forms. Figure 1 shows the 95.45 confidence interval (corresponding to two standard deviations) of the peak error for different noise levels . The results reveal good performance for all methods, all with less than 4 error even at 20 noise level. Although the integrated and the projected‐stress Poisson estimators display a similar variability, the error curve of the PPE is somewhat shifted downwards. This is the case even for , that is, the reason for the better performance of the PPE stems from its smaller discretization error. One possible explanation might be given in a very recent study by Araya et al. : their theoretical estimates describe a suboptimal asymptotic convergence rate for the PPE—although it is still unclear whether their analysis is sharp. For the other estimators, there are no significant changes whether using the Laplacian or rotational description of the viscous terms. The main reason for the small error levels in this example is the fact that the convective term has been “switched off”. When convection is included, the estimated pressure becomes visibly more sensitive to noise, as shown in Figure 2. This indicates that the quadratic nature of the convective term may be more harmful to the estimators' accuracy than the higher‐order derivatives of the viscous term—even for a flow regime with pronounced viscous effects such as this one. This difference will be even more expressive when considering the divergence form of the convective term, as shown in the next test case.

FIGURE 1

Womersley flow benchmark: noise sensitivity (95.45 confidence interval from 30 realizations) considering different versions of the Poisson, Stokes and virtual work‐energy relative pressure estimators

FIGURE 2

Womersley flow benchmark: relative pressure drop (20 noise level, 95.45 confidence interval) when including (lower row) or omitting (upper row) the convective term. The wider bands obtained when considering convection reveal a higher noise sensitivity stemming from convective terms than from viscous terms

Two‐dimensional radial flow

The next test case considers a radial flow that solves both the Navier–Stokes and the Euler equations. This solution describes a pulsatile flow accelerating through a narrowing channel. The velocity and mean pressure drop between inlet and outlet are: where and are the inlet and outlet radii, respectively, and is a given pulsatility, which is chosen here as: The parameters set for this example are  m/s,  rad/s,  kg/m,  mm, and , with an arc of rad. The discretized geometry is shown in Figure 3 (the flow happens from right to left), and the temporal resolution is again  s. Now, the viscous term is dropped in the estimators so as to compare them solely with respect to the treatment of the convective term. The variants considered are the WERP, IMRP, STE, STE, PPE and PPE (cf. Table 1). The differences are now more expressive than in the previous case study: each of the methods is visibly more sensitive to noise when the divergence form of convection is used. Also notice that the errors here are at least one order of magnitude higher than in the previous example, once more indicating the higher noise sensitivity caused by convective (in comparison to viscous) terms. It is important to remark that, in practice, 4D flow data are acquired spatially in a Cartesian‐like manner, so that boundaries will in general not be resolved as in Figure 3; a more realistic setup will be considered in the next test case (Figure 4).

FIGURE 3

Radial flow benchmark: geometry and discretization used for the noise sensitivity study (flow direction: right to left)

FIGURE 4

Radial flow benchmark: noise sensitivity ( confidence interval from 30 realizations) considering different versions of the Poisson, Stokes and virtual work‐energy relative pressure estimators. The conservation form of the convective effects (IMRP, STE, PPE) is more sensitive to noise when compared to the standard (convective) form, for all three classes of estimators

In‐vitro flow in giant intracranial aneurysm phantom

The next test case considers a more realistic scenario using real 4D flow MRI measurements. For an in‐depth study on the flow mechanism in intracranial aneurysms, Amili et al. constructed a 3D‐printed replica of a giant, pre‐treatment aneurysm in an internal carotid artery (ICA). The geometry was scaled (enlarged) by a factor of 2.0 to enable more accurate measurements, and the remaining parameters were adjusted to have an in‐vitro experiment matching the physiological Reynolds and Womersley numbers. The experiment had a velocity encoding of 0.8 m/s along all directions, a flip angle and a spatial resolution of  mm. The pulsation period was  s, with 16 acquisitions per cycle ( s). The fluid properties were measured as  kg/m and  mPa s, and further details on the experimental setup can be found in the original article.

The first step for the pressure computation is to extract a finite element mesh from the data. For this, a tensor‐product hexahedral grid is first constructed directly from the voxelated data. Then, the elements corresponding to voxels outside the aneurysm are deleted (the published data contain a mask assigning a “number/not‐a‐number” flag to each data point). Finally, the mesh is truncated at cm to eliminate the artificial prolongations of the artery used in the experiment (cf. Figure 1 therein ). No erosion or deletion of wall voxels was performed. The final hexahedral mesh, depicted in Figure 5A, contains approximately 2 × 105 elements and nodes. A subsampled version with half of the original spatial resolution is also considered for the last test, see Figure 5B. The relative pressure fields computed at systole by the pressure Poisson and Stokes estimators are illustrated in Figure 6, which shows—as usual—that the PPE estimates a lower relative pressure.

FIGURE 5

Hexahedral meshes constructed from the in‐vitro data by Amili et al.

FIGURE 6

Relative systolic pressure distribution computed from 4D flow measurements in an in‐vitro ICA aneurysm using the PPE and STE. As often observed, the Poisson estimator computes lower relative pressure magnitudes

Figure 7 depicts the mean pressure drop between the inlet and outlet planes. The Stokes and virtual work‐energy estimators show excellent mutual agreement, with less than difference in the peak systolic value. The same goes for the IMRP and STE, both of which use the same integrated version of the convective term. Although there are no reference measurements here, Bertoglio et al. also noted that the (integrated) divergence form of convection increases the estimated pressure for the STE. The results in Figure 7 confirm that this is also the case for the other two families of estimators. The WERP and WERP estimates are virtually indistinguishable, which confirms that, as often assumed in the literature, , , , the viscous boundary term is negligible also for this test case. There is also hardly any difference between the projected‐stress and integrated Poisson estimators, which is most likely due to the rather high resolution of the data (the PPE typically has a higher truncation error, see Figure 1 [right] and also Ref. 46).

FIGURE 7

Pressure drop computed from 4D flow measurements in an in‐vitro ICA aneurysm. The Stokes and WERP estimators agree very well with each other, and so do the IMRP and STE (since both use the same integrated form of the convective term). There is virtually no difference between the WERP and WERP estimates, which confirms that the surface integral is negligible for this flow. Also between PPE and PPE there is hardly any difference, whereas using the divergence form of convection (PPE) increases the pressure drop, as in the other two families of estimators

Next, as done by Nolte et al., noise is artificially added to the data in order to compare the estimators with respect to the different treatment of the convective term. The results, depicted in Figure 8, show rather small variability even for the high noise level applied (). This is due to the relatively high resolution of this dataset, which considers an enlarged aneurysm phantom (error bands for in‐vivo measurements are typically higher). For this reason, no clear winner in terms of noise sensitivity can be picked from the results in Figure 8, with all variants displaying a comparable variability. Then, to yield a more realistic comparison, a spatial subsampling of the original data is applied, with half of the original resolution (cf. Figure 5B). The results, depicted in Figure 9, are now more conclusive. At the pressure peak, the standard deviation of the STE is 62 larger than that of the STE, and the difference is even larger (106) between the IMRP and WERP estimators. For the PPE, though, the divergence form of convection only led to a 16 increase in the standard deviation. These results confirm the trend observed in the in‐silico test case: the divergence form of the convective term is indeed more sensitive to noisy velocity data.

FIGURE 8

Pressure drop ( confidence interval from 30 realizations) computed from 4D flow measurements in an in‐vitro ICA aneurysm, with artificially added noise (). The divergence variants of the convective term (right column) estimate a higher pressure drop in all cases, but there is no “clear winner” with respect to noise sensitivity

FIGURE 9

Pressure drop ( confidence interval from 30 realizations) computed from 4D flow measurements in an in‐vitro ICA aneurysm, with artificially added noise () and spatial subsampling, cf. Figure 5B. The divergence variants of the convective term (right column) estimate a higher pressure drop in all cases, but are more sensitive to noise

DISCUSSION

Main findings

This study compared three families of relative pressure estimators (Poisson, Stokes and virtual work‐energy balance) with different ways of treating viscous and convective terms. Classical treatments as well as two new ones were investigated. The in‐silico tests indicate that, for all three classes considered here, the integrated forms (PPE, STE, and IMRP) yield larger errors for noisy data. Another key finding here is that the noise sensitivity of the convective term can be considerably higher than that of the viscous term (despite the latter's second‐order derivatives) even for cerebrovascular flow regimes with increased viscous contribution. Moreover, the investigation with in‐vitro 4D MRI data confirmed the well‐known underestimation bias of pressure Poisson estimators (in different variants). Notice, though, that new ideas and techniques have recently emerged to increase the PPE's accuracy. For instance, Araya et al. noted that increasing the polynomial order of the pressure interpolant can mitigate the underestimation of peak systolic pressures. A similar effect can be achieved by removing the outer layer of voxels, as recently reported by Nolte et al. A more accurate computation of the viscous term can also be achieved by reconstructing the velocity field near the walls.

Treatment of viscous terms

The (herein introduced) projected‐stress pressure Poisson estimator (PPE) improves the accuracy of the viscous pressure drop, although this improvement is likely to be overshadowed by convection‐related effects in clinical applications. However, the viscous component per se can be a potential biomarker for aortic diseases, in which case a more accurate computation of viscous effects becomes relevant. In particular, the near‐wall velocity reconstruction technique proposed by Donati et al. can be combined with the present viscous modifications to better capture viscous pressure drops. In such a context, the boundary‐integral‐free modification proposed herein for the virtual estimators may prove advantageous. Although whether it could have an impact in clinical practice is unclear at the moment, the WERP estimator has no added cost or complexity in comparison to the classical versions, and has the advantage of not requiring prior phenomenological knowledge on the flow regime by the user, since the boundary integral vanishes exactly for any flow.

Treatment of convection terms

The convective effects in all estimators can be treated in either convective (standard) or divergence/conservation form. Using the latter in Stokes and work‐energy estimators enables integration by parts, thus eliminating derivatives. Different test cases reveal increased noise sensitivity when shifting from the convective to the conservation form—in all three families of methods considered. On the other hand, the study with in‐vitro data showed a larger pressure drop estimated by the conservation variants. Although there are no reference measurements here to judge such increase in estimated pressure, this behaviour is in line with the findings of Bertoglio et al., who reported that the integrated variants can mitigate the pressure underestimation caused by low‐resolution data.

Limitations

The present work is application‐driven, yet still fundamental and, up to a certain extent, also theoretical. The test cases reveal insightful information and trends, but the absence of reference in‐vivo data limits the reach of some of the conclusions drawn here. Moreover, the improved computation of viscous terms does not solve the dependency of relative pressure field estimators (PPE, STE) on the defined geometrical domain. , In fact, viscous effects are most significant near vessel walls, which are often not (accurately) included in the final computational domain. , Therefore, the methods presented here could and should be combined with techniques aimed at improving near‐wall velocity accuracy.

CONCLUDING REMARKS AND FUTURE DEVELOPMENTS

The purpose of this work has been twofold: from one side, to compare different treatments of viscous and convective effects in relative pressure estimation and, from another, to propose modifications in the viscous terms to circumvent surface integrals. The in‐silico and in‐vitro datasets considered here lie in the typical physiological regime of cerebrovascular flows, where viscous forces are more significant than in, for example, aortic applications. In the aorta, it is known that viscous terms are usually negligible if compared to inertial ones, , therefore noise‐related errors due to viscous second‐order derivatives are also typically small (see e.g., Figure 4 in Ref. 1). What is now shown here is that this applies also to flows with lower Reynolds and Womersley numbers. In fact, even in a Womersley flow with no (physical) convective contribution whatsoever, viscous effects can be the least significant source of error. Moreover, using the integrated forms of Poisson, Stokes or virtual work estimators tends to increase noise sensitivity, in spite of their reduced order of differentiation on the velocity data.

As for the modifications introduced here, the in‐silico tests showed that the projected‐stress PPE can improve the accuracy of the viscous pressure drop. It has also the advantage of allowing for local viscosity variations, which could improve the computation of turbulence‐related pressure drops stemming from eddy viscosities. The present rotational WERP, on the other hand, did not provide tangible differences with respect to the standard variant. It can still be seen, however, as more robust in theory, by consistently eliminating (rather than neglecting) viscous surface integrals regardless of the flow regime—without making any assumptions on the relative magnitudes of different virtual work terms. Whether this will be advantageous for flows of clinical interest, however, is not known at the moment.

Two ramifications are planned for the continuation of this study. The first and obvious one is to assess the modified variants presented here in combination with in‐vivo data. The second one involves a more technical aspect. One of the reasons normally attributed to the underestimation bias of pressure Poisson estimators are their increased regularity (continuity) requirements on the pressure field. In this context, on‐going work includes the implementation of an ultra‐weak finite element method for the PPE. In such a framework, integration by parts can be applied to reduce the pressure regularity requirements down to the same ones as in the Stokes estimator, thereby potentially increasing the accuracy of the PPE in low‐resolution grids.