Model Verification and Error Sensitivity of Turbulence-Related Tensor Characteristics in Pulsatile Blood Flow Simulations

Abstract

Model verification, validation, and uncertainty quantification are essential procedures to estimate errors within cardiovascular flow modeling, where acceptable confidence levels are needed for clinical reliability. While more turbulent-like studies are frequently observed within the biofluid community, practical modeling guidelines are scarce. Verification procedures determine the agreement between the conceptual model and its numerical solution by comparing for example, discretization and phase-averaging-related errors of specific output parameters. This computational fluid dynamics (CFD) study presents a comprehensive and practical verification approach for pulsatile turbulent-like blood flow predictions by considering the amplitude and shape of the turbulence-related tensor field using anisotropic invariant mapping. These procedures were demonstrated by investigating the Reynolds stress tensor characteristics in a patient-specific aortic coarctation model, focusing on modeling-related errors associated with the spatiotemporal resolution and phase-averaging sampling size. Findings in this work suggest that attention should also be put on reducing phase-averaging related errors, as these could easily outweigh the errors associated with the spatiotemporal resolution when including too few cardiac cycles. Also, substantially more cycles are likely needed than typically reported for these flow regimes to sufficiently converge the phase-instant tensor characteristics. Here, higher degrees of active fluctuating directions, especially of lower amplitudes, appeared to be the most sensitive turbulence characteristics.

1. Introduction

Flow-phenotypes associated to highly disturbed (chaotic and irregular) hemodynamics play a essential role in the initiation and progression of many cardiovascular disease [1,2,3]. These so-called turbulent-like conditions often prevails at valvular/vascular malformations, but have lately also been found or presupposed under apparent normal physiological flows [4,5]. To-date, the most common modalities to non-invasively estimated these flow conditions are via high-fidelity magnetic resonance imaging (MRI) techniques [6,7] or computational fluid dynamics (CFD) simulation methods [8,9]. Appropriate verification, validation and uncertainty quantification in (image-based) patient-specific cardiovascular numerical modeling are essential steps in order to approach clinical utility [9,10,11,12]. State-of-the-art modeling praxis has been well covered for laminar vascular flows [8,12,13], while general guidelines related to numerical predictions of turbulent-like hemodynamics have received less attention, in spite of the growing number of published turbulence-related CFD studies within the research community. In a series of hemodynamic CFD Challenges [14,15,16], substantial variability in modeling strategies were observed among different research groups worldwide, contributing to widespread result predictions. Similar, large-scale studies in more disturbed turbulent-like flow have still not been initiated. However, the wealth of turbulence-related research studies found today suggests a large gap of consensus among CFD practitioners, with large variability in modeling approaches and verification assessments. Several recent studies have pointed to the importance of having adequate numerical solution strategies (e.g., proper spatial and temporal resolution as well as non-dissipative/diffusive discretization methods) to be able to capture weaker transitional flow regimes as well as high-frequency content in more developed turbulent flows [17,18,19,20,21]; which commonly used CFD methods may fail to predict. Part of these shortcomings may arguably be related to insufficient verification procedures. Model verification makes sure that the numerical model of the particular problem is solved correctly and accurately [22], where recognizable (epistemic) errors related to, for example, discretization and statistical convergence are quantified and evaluated before engaging in the actual CFD analyses. In computational hemodynamics, model verification is often performed using simple, lower-order parameters such as first moments (e.g., mean velocity magnitude), which typically are evaluated at a few spatial locations (points) and/or aggregated over the domain (spatial-averaged) and the cardiac cycle (time-averaged). At the same time, subsequent analyses may involve detailed evaluations of more complex parameters derived from higher moments (e.g., Reynolds stresses). Such verification studies may manifest in misleading results and conclusions.

Periodic pulsatile flows with pronounced (random) cycle-to-cycle flow variations are governed by turbulent-like behavior. A common way to describe these features is by assessing the Reynolds stress tensor, which here describes the momentum flux due to turbulent motion in different directions from a phase-averaged sense. Markers derived from this tensor have, for example, been used for flow-induced blood trauma predictions [23,24,25], turbulence-induce wall effects [26,27,28], and non-invasive assessments of vascular and intraventricular severity [29,30,31,32]. An extensive way to characterize the local properties of such tensors (i.e., second-order symmetric) is to quantify its magnitude (turbulence kinetic energy, TKE) as well as mapping of the anisotropic stress-states [33,34,35]; descriptors which are invariants and therefore well suited for complex geometries associated to cardiovascular flows. The stress-states provide insight into the relative strength of the tensor principal eigenvalues, called ”componentality”, that is, the degree and nature of the turbulence anisotropy, which can be outlined in an invariant anisotropy map (AIM) (Figure 1a). In a recent study, we showcased these tensor characteristics in the turbulent region derived from well-resolved patient-specific CFD simulations [36].

The present study aimed to investigate the Reynolds stress tensor characteristics impact from different CFD modeling strategies conventionally used for verification of pulsatile blood flow predictions, that is, the spatiotemporal resolution and phase-averaging sampling size. These properties were evaluated by assessing the tensor magnitude and anisotropic stress-states in an aortic coarctation (CoA), which was viewed as a representative example of severe patient-specific turbulent flow conditions. With this suggested framework, we hope to pave the way for more complete and robust verification procedures of these flow regimes within the biofluid community, as well as shed some light on the general modeling sensitivity of each turbulence descriptor described herein.

2. Methods

In periodic pulsatile flows a velocity signal $u_i$ in a point $\mathbf{x}(x,y,z)$ can be decomposed into three parts:

$$u_i\left(\mathbf{x},t\right)={\overline{u}}_i\left(\mathbf{x}\right)+{\tilde{u}}_i\left(\mathbf{x},t\right)+u_i^{\prime }\left(\mathbf{x},t\right),$$

(1)

represented by the deterministic mean value ${\overline{u}}_i$ and cycle-related variations ${\tilde{u}}_i$, and the quasi-deterministic turbulence-related fluctuations $u_i^{\prime }$. The six independent phase-averaged correlations of these fluctuations can be used to define the symmetric (3 × 3) Reynolds stress tensor:

$$\begin{array}{cc}\hfill & R_{ij}=\rho \langle u_i^{\prime }{u}_j^{\prime }\rangle =\frac{\rho }{N}\sum _{n=0}^{N-1}\left(u_i-\langle u_i\rangle \right)\left(u_j-\langle u_j\rangle \right),\hfill \\ \hfill & \langle u_i\rangle ={\overline{u}}_i+{\tilde{u}}_i,\hfill \end{array}$$

(2)

where $\rho$ is the fluid density and $[n,N]$ the ensembled cycles range. The anisotropy Reynolds stress tensor is given by:

$$b_{ij}=R_{ij}/2k-{\delta }_{ij}/3=v_{ik}{\mathsf{\Lambda }}_{kl}{v}_{jl},$$

(3)

where $k=R_{kk}/2$ is the TKE. $v_{ij}$ and ${\mathsf{\Lambda }}_{kl}$ are the eigenvector matrix and diagonal eigenvalue matrix, respectively. On canonical form (${\lambda }_1⩾{\lambda }_2⩾{\lambda }_3$), these normalized eigenvalues can be used to construct the barycentric coordinates ($x_B,y_B$) of the AIM:

$$\begin{array}{cc}\hfill & x_B=C_{1C}{x}_{1C}+C_{2C}{x}_{2C}+C_{3C}{x}_{3C},\hfill \\ \hfill & y_B=C_{1C}{y}_{1C}+C_{2C}{y}_{2C}+C_{3C}{y}_{3C},\hfill \\ \hfill & \left\{C_{1C},C_{2C},C_{3C}\right\}=\left\{{\lambda }_1-{\lambda }_2,2({\lambda }_2-{\lambda }_3),3{\lambda }_3+1\right\},\hfill \end{array}$$

(4)

where $C_{iC}$ are weights ($[0,1]$ and $\sum C_{iC}=1$) that measures the closeness to the special limiting states of turbulence anisotropy: one-component (1C), two-component axisymmetric (2C) and three-component isotropic (3C) turbulence. A colorized AIM can be created by associating desirable color triplets to these weights:

$${\left[RGB\right]}^T=C_{1C}\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]+C_{2C}\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]+C_{3C}\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right],$$

(5)

where each realizable turbulence state now corresponds to a specific color, for example, red-green-blue (RGB) that was used in the current study (Figure 1b).

The patient-specific model was reconstructed from MRI data, with valid patient and ethical consents. The ascending aorta inlet pulsatile flow rate were reconstructed from two-dimensional cine phase-contrast MRI measurements assuming a flat (plug-shaped) velocity profile. The aortic arc branching outflows were governed by a square law [37], and the descending aorta outlet assigned with a zero static pressure condition. The fluid viscosity followed a shear rate dependent non-Newtonian relationship [38,39]. For more details the reader is referred to previous studies [27,28,40].

Numerical solutions were computed using the element-based fully coupled (implicit) finite volume solver ANSYS CFX (ANSYS Inc, Canonsburg, PA, USA), using large eddy simulations (LES) and the wall-adaptive local eddy-viscosity (WALE) subgrid model. The discretization of the governing equations was done with the solver’s recommended second-order accurate schemes on high-quality hexahedral cells to minimize dispersion/dissipation errors. ANSYS CFX adopts a median dual mesh strategy (vertex-centered method), where finite-element shape functions are used to approximate the solution gradients at the control volume surfaces. Here, the convective and diffusion terms were computed using (true) tri-linear shape functions (analog to a central difference scheme), whereas a linear-linear interpolation scheme was applied for the pressure gradients. To minimize pressure-velocity decoupling effects, Rhie-Chow interpolation was adopted [41]. The transient term was approximate by a second-order backward Euler scheme using adaptive time-stepping bounded by the Courant–Friedrich–Lewy (CFL) criteria. Within each time-step, iterative solution convergence of each discretized equations was controlled by a domain-based root-mean-square (RMS) residual of 10${}^{-5}$ (with a minimum of 3 iterations) and ensuring a global imbalance target below $1\%$.

It is generally recommended to keep the CFL number below unity (CFL < 1) in every mesh node to ensure appropriate temporal resolution with respect to the local mesh size and fluid velocity. Depending on the flow case, however, this restrict time-step criteria can contribute to substantial over-resolved temporal characteristics in large parts of the computational domain. In our previous studies of the same numerical model [27,40], a more relaxed CFL criteria of maximum CFL < 5 captured near similar cardiac evolution of the TKE field in comparison to the maximum CFL < 1, with and without adaptive time-stepping. In this case, a maximum CFL < 5 resulted in global mean CFL ∼ 0.5. In this study, a adaptive time-stepping profile was prescribed into the solver, derived from the built-in adaptive scheme in the CFD solver. In

Table 1. Temporal statistics of the different time-stepping approaches, including the range of the time-step sizes ($\mathsf{\Delta }t$, in microseconds $\mu s$), number (#) of time-steps and computational costs (in CPU-hours) per cycle.

	3MC, maxCFL < 5	3MC, maxCFL < 5	6MC, maxCFL < 1
	(Adaptive $\mathbf{\Delta }\mathit{t}$)	(Constant $\mathbf{\Delta }\mathit{t}$)	(Adaptive $\mathbf{\Delta }\mathit{t}$)
$\mathsf{\Delta }t$ range, [min, max]	[162, 2000]	162	[20, 771]
# of $\mathsf{\Delta }t$ per cycle	2000	6200	14,000
CPU-hours per cycle	660	2000	8700

, the temporal statistics and run-time costs from the adopted time-stepping schemes are shown, for perspective also including the conventional constant time-step approach.

Sensitivity Analyses

The spatial–temporal sensitivity analysis were evaluated by comparing a model represented by 6 million cells (MC) with maximum CFL < 1 (6MC-CFL1) against a 3 MC with maximum CFL < 5 (3MC-CFL5). This corresponded to an averaged grid refinement factor of roughly 1.3 in the disturbed turbulent flow region, which may be viewed as sufficiently fine according to general recommended procedures for assessing CFD-related discretization errors [42]. Fifty cardiac cycles ($N=50$) were used to reconstruct the turbulence statistics in both cases, corresponding to a computational cost of around 480 k and 36 k (i.e., ∼13 times less) CPU-hours for the 6MC-CFL1 and 3MC-CFL5 case, respectively, using a state-of-the-art supercomputer (National Supercomputer Centre, Linköping, Sweden). This was viewed as an upper limit concerning computational resources.

For the phase-averaged sensitivity analysis, a total of 80 cardiac cycles were acquired using the 3MC-CFL5 case, where different non-overlapping cycle ranges were compared against each other (denoted * vs. **) as well as against the whole data range. For example, 20* vs. 20** compare the phase-averaged results between the first 40 cycles against the last 40 cycles, whereas for example, the notation 40* vs. 80 compares the first 40 cycles against the complete range ($N=80$). Here, five initial cycles were simulated before collecting data in all studies to minimize initialization effects on the results. Data was saved at 100 equally spaced time-steps every cardiac cycle.

The level of agreement between different cases was evaluated, as an example, across several post-stenotic cross-sectional planes, both qualitatively as well as quantitatively, using a point-based root-mean-square deviation (RMSD) measure, which for the componentality weights $C_{iC}$ was defined as:

$$C_{rms}=\sqrt{\frac{1}{3}\sum _{i=1}^3{\left(\mathsf{\Delta }{C}_{iC}\right)}^2},$$

(6)

where $\mathsf{\Delta }$ represents the cell node difference between two cases. The corresponding RMSD of TKE is simply the absolute difference ($k_{rms}=\sqrt{{\left(\mathsf{\Delta }k\right)}^2}=\left|\mathsf{\Delta }k\right|$). Here, $C_{rms}$ provides a convenient percentage sense of the relative deviation of the stress-states within the AIM. To attain an equivalent representation, $k_{rms}$ was also normalized by the third quartile ($Q_3$) TKE value in the related turbulent region (i.e., the median of the upper half of the dataset, or 75th percentile) to put more emphasize against higher TKE intensities. The $Q_3$ spatiotemporal-averaged value across all cross-sections was also used to evaluate the general tendency of each parameter against the phase-averaging sample size.

3. Results

This section will start with a concise description of the general flow characteristics in the post-stenotic region. For more details, the reader is referred to earlier studies [28,36]. Four main flow stages were identified throughout the systolic part of the cardiac cycle. At the first two stages (mid-acceleration to early-deceleration), an eccentric jet was formed, followed by shear-layer destabilization and transition to turbulence. Over the third, early flow deceleration (EFD) stage, quasi-steady post-stenotic flow patterns could be noticed, with evident vortical breakup, and turbulence intensification in the trace of the jet. At the succeeding late flow deceleration stage, the post-stenotic jet began to collapse, followed by turbulence relaminarization. This study focused on the EFD stage, where the turbulent instabilities were most profound. For reference, high TKE intensity is herein referred to values above 60 Pa (i.e., purple and higher in the range of $[0,160]$ Pa colorbar), which corresponds to the around the upper $25\%$ (third quartile) of the post-stenotic TKE values.

3.1. Spatiotemporal Sensitivity Analysis

The phase-instant Reynolds stress characteristics showed a general qualitatively agreement between the two resolutions (Figure 2), however, also with local spots of substantial deviation (TKE > 50 Pa, $C_{rms}$ > 0.25). The turbulence states (Figure 2b) appeared less coherent in comparison to the TKE distribution (Figure 2a). Generally, the more anisotropic regions (1C and near-wall 2C) coexist best between the cases. No clear association could be seen between elevated $C_{rms}$ and $k_{rms}$ regions. The spatially averaged (mean) errors were fairly consistent (∼15%) between the two metrics (Figure 3, Instant, black bars).

From a time-averaged perspective (Figure 4), similar but more coherent overall Reynolds stress characteristics were observed, owing to the near-steady flow patterns over this period. Here, a much stronger overall agreement could be noticed, with a more than a twofold and fivefold reduction of the most extreme RMSD value for the TKE (Figure 4a) and stress-states (Figure 4b), respectively. Collectively, the cross-sectional mean error was reduced to ∼5% for both metrics (Figure 3, EFD, green bars).

3.2. Phase-Averaging Sensitivity Analysis

Verification of the phase-averaging sampling size was performed across several cross-sectional locations downstream the constriction (0D to 3D), from a phase-instant and time-averaged perspective. The overall findings at these locations were similar; therefore, only qualitative results from one cross-section (1D) were presented (Figure 5 and Figure 6), whereas the deviation errors and $Q_3$ trends, were considered as mean values over all cross-sections (Figure 3 and Figure 7).

For the phase-instant results (Figure 5), the elevated ring-like TKE structure surrounding the jet appeared more coherent first at the 10 cycle range, however, not without large local deviations (RMSD TKE $>$ 80 Pa); also evident for 20 cycles. Up to the 20 cycle range, the mean TKE error (Figure 3, Instant, blue and brown bars) were well $>$20% for the independent ranges (* vs. **) and $>$10% against the 80 cycles reference case. Concerning the stress-states, 10 cycles or lower resulted in very poor qualitative resemblance compared to the full range. At 20 cycles, parts of the anisotropic 1C-like could be depicted, but still with overall large RMSD. Comparing 40* vs. 80, the local RMSD and mean errors were clearly reduced and similar ($\sim$8%) between the metrics. Here, however, the two independent 40 cycle ranges did not fully correlate statistically, exhibiting a mean error $>$12% for both TKE and stress-states.

For the time-averaged analysis (Figure 6), the difference in Reynolds stress characteristics was generally much lower than the phase-instant results. In comparison to the 80 cycles, the TKE patterns showed generally good qualitative agreement as low as for 10 cycles, with a mean error ∼7% (Figure 3, EFD, gray and red bars), while a similar level of agreement for the stress-states required 20 cycles or more. To reach a mean error below $5\%$, 20 and 40 cycles were required for the TKE and stress-states, respectively. Here, however, it is also important to acknowledge that surprisingly small RMSD (or strong resemblance) between the lower independent cycle ranges, despite the poor agreement against the reference case. Evident is also the movement from 1C-like turbulence towards 2C state in the near-wall region as the number of cycles were increased in the phase-averaging procedure.

The sensitivity on the third quartile value $Q_3$ (Figure 7), averaged over the EFD phase, pointed out a general underestimation trend of the high TKE and $C_{3C}$ values for the reduced cycle ranges, whereas the $C_{1C}$ metric was generally overestimated. The $C_{2C}$ metric was clearly underestimated in the 5–10 cycle range, while for more cycles, a mixed tendency could be seen, consequently with an unreliable mean $Q_3$ value due to cancellation effects. Compared to the 80 cycles, a near $5\%$ mean error margin was attained at 20 and 40 cycles for the TKE and $C_{1C}$, respectively. At 40 cycles, however, the weight governing the isotropic state ($C_{3C}$) still showed large discrepancies (∼40%) against 80 cycles; findings which could also be noticed in the 40* vs. 80 AIM (Figure 6). At this range the $C_{2C}$ min/max deviation was still considerably large ($[-28\%,+24\%]$).

4. Discussion

In this work, the sensitivity of two important CFD modeling-related solution errors associated with the Reynolds stress characteristics were explored in the post-stenotic region of a patient-specific coarctation model using scale-resolving simulations. The proposed verification method considered both the magnitude and nature of the turbulence anisotropy; flow properties that have previously shown to vary substantially in these flow regimes [28,36]. In summary, these turbulence descriptors appeared most sensitive for the phase-instant compared to time-averaged results. For a low amount of cycles, the overall phase-averaged error clearly trumped the errors related to the spatiotemporal resolution of the computational model. Overall, the turbulence componentality was the more sensitive property to satisfy statistically, especially the more isotropic stress-states of moderate-to-low turbulence intensity.

In computational hemodynamics, great emphasis is usually put on choosing proper discretization schemes, mesh, and time-step resolution to predict turbulent-like flow characteristics [12,13,17,18,19,21], meanwhile phase-averaging aspects have not received much attention. This oversight may partly be explained by the relatively large computational costs associated with patient-specific pulsatile simulations compared to laminar or steady flows, making these types of sensitivity analyses very expensive/unpractical for the average CFD practitioner. Of the studies associated to these flow regimes, the range of cycles used for phase-averaging seems to vary substantially, with for example, reported ranges from ∼5–10 cycles [31,32,43,44,45,46,47,48], ∼20–50 cycles [36,49,50,51,52,53,54] and up to ∼80–100 cycles [28,55]. From a statistical point-of-view, data consisting of even 100 samples seems to be on the borderline to capture the strongly irregular behaviors in turbulent flows. However, due to the periodic nature of these conditionally turbulent flows, the variance of the cycle-to-cycle flow fluctuations are most likely more constrained in comparison to for example, steady flow conditions where more broadband turbulence characteristics can be expected [36,54]. Although, to obtain sufficiently low phase-averaged related errors in these lower cycle ranges presumably involves analyzing crude hemodynamic quantities (e.g., substantially aggregated lower-order moments).

The present study showed that the phase-averaged error induced from using too few cycles could easily outweigh the error associated with spatiotemporal resolution. For example, an error margin $<$10% could not be satisfied at 40 cycles for either metrics (Figure 3, Instant), that is, also without including the spatiotemporal error. Here, the phase-instant errors between the spatiotemporal resolutions (acquired over 50 cycles) and the 40* vs. 40** independent cycle ranges were comparable (∼14–16%), which may suggest that the former is governed mainly by the latter errors. These phase-relate errors could almost be halved by doubling the cycle range (40* vs. 80). From a qualitative sense, however, the phase-instant TKE patterns appeared to agree better than the anisotropic states (Figure 5), which still showed relatively fragmented (less coherent) stress-states at these upper cycle ranges. These non-coherent bulk regions are likely associated with the ”slow” phase-averaging convergence of the underestimated isotropic-dominant stress-states, clearly evident for the time-averaged results (Figure 7). In contrast, the more anisotropic states (with less active fluctuation directions) close to the wall and in the jet region appeared to be more forgiving. While similar numbers of assessed cycles have shown to have a small impact on individual Reynolds stress components in controlled pulsatile pipe flows [51] (∼3%, 30 vs. 50 cycles), the findings in this study suggest that at least a twofold more cycles are needed to attained sufficiently small error ($<$ 5%) of the phase-instant tensor magnitude and componentality characteristics.

When time-averaging the results over the most turbulence-prone phase, similar but more coherent TKE and stress-state patterns were observed in both sensitivity studies (Figure 4 and Figure 6). These results may be explained by the quasi-steady jet-like flow over this period, which in effect adds much more samples into the phase-averaging statistics (in this case, 20-times more compared to the phase-instant data). Here, a sufficient phase-averaged error margin ($<$ 5%) for TKE and stress-states could be reached at 20 and 40 cycles, respectively (Figure 3, EFD). The error associated with the spatiotemporal resolution was reduced by a near factor of three (to ∼5%) for both metrics, which may be closer to the error expected in the phase-instant results if much more cycles were to be included. However, it is important to note that the period of these quasi-steady flow conditions may be very case sensitive. Therefore, it is vital to outline the degree of these conditions to assess a proper time-averaging window. Generally speaking, highly aggregated parameters should also be used with causing, as these methods have been shown to be less susceptible to different solution strategies [18]. That said, an interesting investigation would be to incorporate more data into the time-averaging procedure by saving at a higher rate within the cardiac cycle, which currently is several orders of magnitude coarser than the time-step size used in the simulations. Besides data storage challenges, this framework might lead to a noteworthy reduction in required cycles for the phase-averaging verification procedure.

The tendency of the metrics $Q_3$ values (Figure 7) showed a clear $C_{1C}$ overestimation for the lower cycle ranges compared to the 80 cycles reference case, while TKE, $C_{2C}$, and $C_{3C}$, on the other hand, were substantially underestimated; especially the $C_{3C}$ weight. These trends can qualitatively also be appreciated in the contour plots (Figure 6), where 1C-dominant characteristics occupied a significant portion of the cross-sections for the lower cycle ranges. These findings are reasonable as the random characteristics governing the $C_{3C}$ isotropic state require more samples to statistically satisfy, as opposed to stress-states with less active components ($C_{1C}$ and $C_{2C}$). Indeed, even at the 40 cycle range, $C_{3C}$ was still far from converged, with a deviation ∼40%. However, considering the low deviations for TKE-based $Q_3$ value already at 20 cycles (∼5%), this lack of convergence is mostly associated with moderate-to-low (first and second quartile) TKE intensities, which qualitatively can be recognized in the stress-state contour plots (Figure 5 and Figure 6). Contrarily, Reynolds stress characteristics manifested by high TKE values and degree of anisotropy were easier to converge from a phase-averaging sense. Likewise, in the near-wall anisotropic region, that is, stress-states close to the two-component limit (Figure 1), the two lower cycle ranges showed a clear tendency to overestimate the 1C characteristics (Figure 6), while approaching more axisymmetric $C_{2C}$ characteristics as more cycles were included. These observations may partly also explain the $C_{1C}$ and $C_{2C}$ general $Q_3$ trends in these lower ranges (Figure 7). A final remark was also the small error between each independent time-averaged cycle range (Figure 3, EFD, * vs. **), which suggests that this verification strategy cannot be used alone, without also comparing against a much more resolved case. Indeed, it is first, when both these verification strategies satisfy the desirable error-criteria, that the solutions should be considered converged, which inherently is fulfilled for a large number of independent samples, but, as shown here, clearly not for too low amounts on specific flow descriptors. In this case, a possible way around these issues may be to perform the verification analyses merely on the most sensitive parameter (i.e., $C_{3C}$ weight).

This work contains several assumptions related to the patient-specific numerical modeling and analyses restrictions, none of which are believed to change the generality of the findings in this study. For detailed insight into specific CFD method assumptions, the reader is referred to our previous studies [27,28,36]. To showcase the verification framework in this study, only one flow case was presented, while substantially more patient-specific turbulent flows are needed to draw more general conclusions. Adding more models were also unpractical from a computational cost standpoint. However, we are still confident that similar error-trends, as seen in these sensitivity analyses, are expected in other pulsatile turbulent-like flows of similar characteristics. Under weaker (less developed), transitional-like conditions, more turbulence anisotropy (or less isotropization) is expected [36,56], which according to the findings herein may require less amount of cycles to statistically converge compared to more well-developed turbulence fields. The true phase-instant error, exclusively associated with the different spatiotemporal resolutions, could not be estimated in this study due to the limited number of used cycles ($N=50$). Indeed, according to the findings herein, this error appeared to be mostly dominated by the phase-averaging procedure, suggesting the spatiotemporal error to be substantially lower. It may also be relevant to investigate the CFL criteria (temporal resolution) impact on these tensor characteristics. While our previous studies only showed minor difference on the overall volume-integrated TKE field between a maximum CFL < 5 and CFL < 1 condition using the same mesh [27,40], the local tensor characteristics might be more sensitive, especially the phase-instant characteristics. It is also important to emphasize that quantitative point-wise deviation measures can also be misleading, where a slight shift in, for example, the jet position (angle, penetration depth, etc.) between two flow cases inherently manifests as large local errors, meanwhile qualitatively the overall flow characteristics may be similar. In this study, the relative modeling errors were, as an example, computed based on RMSD measures on the tensor descriptors and relative to their $Q_3$ values. Of course there are are several other ways to weight these tensor-related errors that can be more suitable depending on the nature of the study. The focus in this work was to showcase new and extensive ways to perform modeling sensitivity analyses in patient-specific turbulence hemodynamics, while also provide some insight into the relative error sensitivity of these descriptors. The errors associated to the turbulence-related tensor properties described herein could also be as estimated by more well-accepted procedures within the CFD community [42], where for example the apparent order of accuracy in the CFD method could be quantified and utilized to estimate a extrapolated relative errors as well as an “convergence index” associated with both the discretization and phase-averaging errors. Interesting continuing work could also be to consider the impact from different pulsatile waveforms (e.g., Reynolds number and Womersley number), high-order numerical schemes, considering much more cycles, a larger sequence of different mesh resolutions, and higher data sampling rate, as well as uncertainty quantification, evaluation of other tensors characteristics (e.g., dissipation tensor), including turbulence structural and orientation information [34,57].

5. Conclusions

This study has established a comprehensive and at the same time practical verification framework for patient-specific CFD predictions of turbulence-related tensor fields, from which various surrogate markers often are extracted. Based on the findings in this work, attention should also be focused on reducing the phase-averaged-related error, as it could easily outweigh the error associated with the CFD model spatiotemporal resolution if too few cardiac cycles are sampled. Concerning the Reynolds stress characteristics, the findings in this study suggest that significantly more cycles, than typically reported for these flow regimes, are needed to statistically convergence the turbulence magnitude and shape of the stress-states (componentality). Here, the phase-instant results were the most sensitive descriptors, particularly moderate-to-low turbulence intensities with more active fluctuating directions, for example, the isotropic dominant states.