2.1 Streamlines through a control volume

Let us consider any flow field $\boldsymbol{v}(\boldsymbol{x},t)$ defined in a spatio-temporal domain ${\mathcal{D}}\times {\mathcal{T}}$ , representing the local fluid speed at any point $\boldsymbol{x}\in {\mathcal{D}}$ at any time $t\in {\mathcal{T}}$ and with a fluid density $\unicode[STIX]{x1D70C}(\boldsymbol{x},t)$ . A streamline at a time $t$ is a line $\unicode[STIX]{x1D6FE}$ whose tangent vector is parallel to the flow field at that time at all points along the line. Parameterized with a pseudo-time $\unicode[STIX]{x1D70F}$ ,

where we consider the interval, $I\subset \mathbb{R}$ , that maximally extends the streamline in ${\mathcal{D}}$ . We will denote by $\unicode[STIX]{x1D6FE}$ (no boldface) the streamline trajectory regardless of its parameterization. The congruence of streamlines, $\unicode[STIX]{x1D6E4}(t)$ , is the set of all streamlines in ${\mathcal{D}}$ at any given instant $t$ . For simplicity, let us fix a time $t$ , and assume it in the notation, as for $\boldsymbol{v}(\boldsymbol{x}):=\boldsymbol{v}(\boldsymbol{x},t)$ and $\unicode[STIX]{x1D6E4}:=\unicode[STIX]{x1D6E4}(t)$ .

Let us consider a control volume ${\mathcal{D}}_{{\mathcal{S}}}\subset {\mathcal{D}}$ enclosed by a control surface ${\mathcal{S}}$ , with a strictly positive scalar field, $\unicode[STIX]{x1D6FC}(\boldsymbol{x})>0~\forall x\in {\mathcal{D}}$ , defining a current $\boldsymbol{j}(\boldsymbol{x})=\unicode[STIX]{x1D6FC}(\boldsymbol{x})\,\boldsymbol{v}(\boldsymbol{x})$ .

In particular, we will consider the volume flow rate, $\unicode[STIX]{x1D711}_{\boldsymbol{v}}^{{\mathcal{S}}}(\boldsymbol{x})$ , corresponding to $\unicode[STIX]{x1D6FC}(\boldsymbol{x})=1$ , and the mass flow rate generated by $\boldsymbol{j}(\boldsymbol{x})=\unicode[STIX]{x1D70C}(\boldsymbol{x})\,\boldsymbol{v}(\boldsymbol{x})$ .

The flux density will be negative at entry points and positive at exit points. Thus, the control surface ${\mathcal{S}}$ can be accordingly split into two regions:

There is indeed a third region, $\overline{{\mathcal{S}}}_{0}:=\{\boldsymbol{x}\in {\mathcal{S}}\,|\,\unicode[STIX]{x1D711}_{\boldsymbol{v}}^{{\mathcal{S}}}(\boldsymbol{x})=0\}$ , where the streamlines are tangent to the surface. In particular, if ${\mathcal{S}}$ includes rigid walls delimiting the fluid, they will be included in $\overline{{\mathcal{S}}}_{0}$ .

Figure 2.

Example of recirculation at the entrance of a control volume ${\mathcal{D}}_{{\mathcal{S}}}$ , represented in red (a). One streamline is shown with recirculation number $r_{{\mathcal{S}}}(\unicode[STIX]{x1D6FE})=2$ . The two inward points are represented in white and the two outward points in black (b). The entrance is divided into two surfaces, $\overline{{\mathcal{S}}}_{In}$ and $\overline{{\mathcal{S}}}_{Out}$ , according to the flux sign.

A streamline $\unicode[STIX]{x1D6FE}\in \unicode[STIX]{x1D6E4}$ can miss ${\mathcal{D}}_{{\mathcal{S}}}$ , enter once or repeatedly enter and exit ${\mathcal{D}}_{{\mathcal{S}}}$ . Let us denote the number of times that the streamline enters and exits ${\mathcal{D}}_{{\mathcal{S}}}$ by, respectively,

The $r_{{\mathcal{S}}}^{In}(\unicode[STIX]{x1D6FE})$ points in the intersection $\unicode[STIX]{x1D6FE}\cap \overline{{\mathcal{S}}}_{In}$ can be sorted according to the crossing order: $\boldsymbol{x}_{In,i}(\unicode[STIX]{x1D6FE})$ , for $i=1,2,\ldots ,r_{{\mathcal{S}}}^{In}(\unicode[STIX]{x1D6FE})$ . Analogously for the $r_{{\mathcal{S}}}^{Out}(\unicode[STIX]{x1D6FE})$ points in $\unicode[STIX]{x1D6FE}\cap \overline{{\mathcal{S}}}_{Out}$ .

Assuming that there is no sink or source inside the control volume, each streamline will in general exit as many times as it enters. Strictly speaking, some exceptional cases can appear, involving equilibrium points or cycles. However, under reasonable assumptions, they would represent a zero-measure set of streamlines, and will be discarded in the following definition.

Figure 2 illustrates these concepts.

This definition guarantees that there exist bijective mappings between inlet, outlet and the congruence of streamlines $\unicode[STIX]{x1D6E4}_{{\mathcal{S}}}$ . Notice that, although the union of the streamlines fills the three-dimensional domain, the congruence of streamlines is a two-dimensional set. Since $\boldsymbol{x}_{In}:\unicode[STIX]{x1D6E4}_{{\mathcal{S}}}\rightarrow {\mathcal{S}}_{In}$ is a bijection, $\boldsymbol{x}_{In}(\unicode[STIX]{x1D6FE})$ can be used as a parameterization of $\unicode[STIX]{x1D6FE}$ . Analogously, $\boldsymbol{x}_{Out}(\unicode[STIX]{x1D6FE})$ gives an alternative parameterization.

Figure 3.

Three examples of streamlines bundles generated by their crossing at the outlet (top) defined as a circle around a point with three different radii. The crossing at the inlet (bottom) is very elongated and the complexity of the obtained pattern depends on the scale considered.

2.2 Probability distribution of streamlines

The computation of any statistical quantity from a streamline distribution requires the definition of a probability density of streamlines. A natural selection is to make the probability density proportional to the local flux, thereby considering denser presence of streamlines where the flow rate is larger.

Any bundle of streamlines, $\unicode[STIX]{x1D6F6}\subset \unicode[STIX]{x1D6E4}_{{\mathcal{S}}}$ , will define a subregion of the inlet by its crossing, ${\mathcal{A}}_{In}(\unicode[STIX]{x1D6F6})=\unicode[STIX]{x1D6F6}\cap {\mathcal{S}}_{In}$ , and a subregion of the outlet, ${\mathcal{A}}_{Out}(\unicode[STIX]{x1D6F6})=\unicode[STIX]{x1D6F6}\cap {\mathcal{S}}_{Out}$ (see figure 3). Accordingly, the probability of a streamline bundle $\unicode[STIX]{x1D6F6}\subset \unicode[STIX]{x1D6E4}_{{\mathcal{S}}}$ will be

Analogously, $p_{Out}(\boldsymbol{x})$ and $P_{Out}(\unicode[STIX]{x1D6F6})$ are defined by their crossing at the outlet. Indeed, for any surface intersecting once every streamline in $\unicode[STIX]{x1D6E4}_{{\mathcal{S}}}$ , we could define the corresponding probability density.

Strictly speaking, we are considering here volumetric probabilities instead of probability densities, in the sense discussed in Tarantola (Reference Tarantola2005) and Tarantola & Mosegaard (Reference Tarantola and Mosegaard2007), since they are defined with respect to the surface element, $\text{d}S$ , and independently of any coordinate system.

Observe that for a homogeneous incompressible flow, the probability distribution generated by the mass flux coincides with the probability distribution generated by the volume flux.

In these cases, the natural probability density is a property of the streamlines, independently of the cross-section considered. We will symbolically denote it as $p(\unicode[STIX]{x1D6FE})$ . Note that the naive homogenous distribution would be cross-section dependent. In general, it would be inhomogeneous for any other cross-section.

2.3 Source of non-invariance of the natural distribution of streamlines

Whenever the current $\boldsymbol{j}=\unicode[STIX]{x1D6FC}\boldsymbol{v}$ is not divergence free, $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{j}\neq 0$ , the corresponding flux cannot be guaranteed to be conserved, and the probability distribution of streamlines will not be invariant. In our experiments, we will consider homogeneous incompressible flows, which should satisfy $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{j}=0$ . However, when a flow is obtained by computer flow simulation, numerical inaccuracies might arise, which justifies studying their influence in the probability distribution invariance.

For divergence-free currents, $\unicode[STIX]{x1D709}(\unicode[STIX]{x1D6FE})=0$ , so that the flux is conserved. Thus, the factor exponent $\unicode[STIX]{x1D709}(\unicode[STIX]{x1D6FE})$ is the source of non-invariance of the corresponding natural distribution of streamlines. For incompressible flows ( $\boldsymbol{j}=\boldsymbol{v}$ ), the integral simplifies to

2.4 Natural probability distribution of world-lines

For unsteady flows, particle trajectories (pathlines) are generally different from streamlines. Pathlines are normally understood as the trajectories in the space, considering the temporal dimension only as the pathline parameter,

However, particle trajectories can be understood as streamlines in the 4-D Galilean space–time, called world-lines. Let us consider the 4-D Cartesian coordinates and the 4-D velocity field

where Greek indices range from 0 to 3. This notation is partially borrowed from special relativity physics, but it is compatible with Galilean relativity and Newtonian mechanics.

World-lines represent the particle trajectories in the space–time and are computed as the streamlines of this 4-D flow:

which, together with the initial condition $\unicode[STIX]{x1D6FE}^{0}(\unicode[STIX]{x1D70F}_{0})=\unicode[STIX]{x1D70F}_{0}$ , is equivalent to the pathline equation

The 4-D nabla operator is denoted

The mass flow current is then $j^{\unicode[STIX]{x1D707}}:=\unicode[STIX]{x1D70C}v^{\unicode[STIX]{x1D707}}$ , for which the continuity equation is simply written in terms of the 4-D divergence:

Hence, all definitions and properties described above for streamlines of incompressible flows are valid for world-lines of any (mass preserving) flow. Thus the natural distribution of world-lines appears as follows. Considering a 4-D control hyper-volume, its boundary would be a 3-D hyper-surface, from which the corresponding inlet and outlet hyper-surfaces can be defined according to the flux density derived from $j^{\unicode[STIX]{x1D707}}$ . This flux defines then a probability density invariant to the cross-sectional hyper-surface.

3.1 Scale-dependent probability model

A typical method to estimate a probability density from a sample is based on binning the variables of interest and extracting the frequencies in each bin (histogram). If the bin size is small compared with the spatial frequency of change of the considered probability density, this discretization provides a valid approximation. However, for the joint probability density, $p(\boldsymbol{x}_{In},\boldsymbol{x}_{Out})$ , defined with infinite precision, any binning would be, by definition, large. Thus, it would impose a coarse-grained scale of observation. This observation scale is, indeed, what we need to introduce in the probability model. Unfortunately, this binning method has sometimes led to confusion between two different concepts: observation scale and numerical estimation. Thus, for clarity, we will explicitly decouple both.

We will introduce the scale by considering that each streamline does not correspond to an exact point at the inlet or outlet, but to a distribution around that point, which we will call point spread function (PSF), noticing its analogy with imaging systems. The PSF is given by independent conditional probability densities at the inlet, $p_{In,s}(\boldsymbol{x}\mid \unicode[STIX]{x1D6FE})$ , and outlet, $p_{Out,s}(\boldsymbol{y}\mid \unicode[STIX]{x1D6FE})$ , providing the joint probability density

where $\boldsymbol{x}\in {\mathcal{S}}_{In}$ and $\boldsymbol{y}\in {\mathcal{S}}_{Out}$ . The index $s$ identifies the scale considered. As introduced in § 2.2, $p(\unicode[STIX]{x1D6FE})$ denotes the natural probability density of streamlines. The integral can be computed either at the inlet or the outlet, parameterizing the streamlines by the corresponding crossing points.

From this joint probability density, we obtain the MI relative to the PSF at scale $s$ :

In contrast with the binning approach, here the scale is explicitly stated as part of the exact definition of the MI, and decoupled from the discretization needed for its numerical estimation.

It is possible to define PSFs equivalent to the binning method, using a discretization of the inlet and outlet. However, in that case, the definition of the joint probability density and its MI will be dependent, not only on the scale but also on the particular partition considered. Besides, the definition of a discretization is not straightforward for inlets or outlets of general irregular shape.

We will consider as PSF, Gaussian kernels with isotropic standard deviation $\unicode[STIX]{x1D70E}$ :

where the normalization constant depends on the inlet geometry. The standard deviation, $\unicode[STIX]{x1D70E}=sR$ , is defined as proportional to the scale, $s$ , and the inlet area-equivalent radius, $R=\sqrt{\text{Area}({\mathcal{S}}_{In})/\unicode[STIX]{x03C0}}$ . For the outlet, $p_{Out,s}(\boldsymbol{y}\mid \unicode[STIX]{x1D6FE})$ is analogously defined. This PSF can be applied independently of the geometry of the control volume, is homogeneous and isotropic by construction and produces continuous and smooth probability densities.

3.2 Interlacing complexity index (ICI)

Several versions of normalized MI for discrete variables exist (Cahill Reference Cahill2010; Vinh, Epps & Bailey Reference Vinh, Epps and Bailey2010; McDaid, Greene & Hurley Reference McDaid, Greene and Hurley2011). They run along the finite range between 0 (if and only if the two variables are independent) and 1 (if and only if there exists a bijection between them). The normalization facilitates interpretation and comparison across different conditions, where the (unnormalized) MI might have different ranges. Two of the normalized MIs are similarity metrics (Chen, Ma & Zhang Reference Chen, Ma and Zhang2009),

and their complements,

are normalized distance metrics (Vinh et al. Reference Vinh, Epps and Bailey2010) (satisfying positive definiteness, symmetry and triangle inequality). Among them, $\mathit{NMI}_{max}$ seems to provide the most linear and intuitive behaviour (McDaid et al. Reference McDaid, Greene and Hurley2011).

Unfortunately, none of the normalized MI versions is directly applicable to continuous variables, since the normalization factor involves either the entropy of each variable or the joint entropy, which are inconsistent for continuous variables. For discrete variables the entropy can be expressed as self-information, $H(\boldsymbol{X})=I(\boldsymbol{X},\boldsymbol{X})$ , but it diverges for continuous variables. However, analogous to (3.2), we can introduce a self-probability density dependent on the observation scale:

and the corresponding inlet self-information, $I_{s}(\boldsymbol{X},\boldsymbol{X})$ , analogously to (3.3), but integrating both variables in ${\mathcal{S}}_{In}$ . In case the scale was introduced by a discretization, it can be proved that $I_{s}(\boldsymbol{X},\boldsymbol{X})$ coincides with the entropy of the corresponding discrete variable. Thus, this allows us to introduce a generalization of the normalization involved in $\mathit{NMI}_{max}$ and $d_{max}$ .

Note that the ICI is not directly the normalized distance $d_{max}$ , but a generalization of it. Thus, we have not proved that it satisfies all its properties. However, the results of our experiments are compatible with the triangle inequality and the expected range, $\mathit{ICI}_{s}\in [0,1]$ . In particular, $\mathit{ICI}_{s}=0$ for the case of a steady flow in a straight cylindrical tube, where the streamlines are straight and parallel. ICI increases with the complexity of the flow, approaching $\mathit{ICI}_{s}=1$ when the positions at the outlet are highly unpredictable from the positions at the inlet.

Figures 7, 9, 11 and 15 show the ICI obtained for different types of geometries and flows, which are described in the following sections. Each of them is presented together with a corresponding illustration of the complexity introduced by the interlacing of the streamlines when propagating from the inlet to the outlet. This provides an intuitive understanding of the meaning of the proposed ICI as a flow complexity quantification.

Observe that all the elements involved in the definition of the ICI are valid for an inlet and outlet of any geometry and topology. In particular, inlet and outlet may consist of multiple, separate regions. which would be analogous to a multiple-inputs multiple-outputs (MIMO) communication system.

4.1 Flow cases

In §§ 5–8, we present several numerical experiments involving five different flow cases. Their characteristics are summarized in table 1. They illustrate a diversity of geometries and flow conditions, including steady and unsteady flows, and different inlet and outlet geometries and topologies. In all cases, except for case (v), rigid impermeable walls and incompressible Newtonian flow is assumed.

Depending on the case, the flow fields and streamlines have been obtained analytically or by numerical simulation and integration. For the simplest case (i), flow field, streamlines and MI have been computed both analytically and numerically. This have been used in § 5 to estimate the accuracy of the algorithms and to test the correctness of their implementation. Cases (ii) and (iii) represent idealized geometries and steady flows, and their corresponding ICI under different conditions is investigated in § 6. In contrast, case (iv) includes anatomically realistic vasculatures and unsteady flows, which are presented and investigated in § 7. A subset of them has been also used in § 5 to evaluate the algorithm accuracy in more complex cases. Case (v) is analytically defined and represents a family of unsteady compressible 2-D periodic cellular flows. In § 8, we analyse the obtained $\mathit{ICI}_{s}$ for different parameter values and test its limit for infinitesimal scales.

Table 1.

Summary of the characteristics of the 4 flow cases considered for the numerical experiments.

4.2 Flow simulations

From the surface triangular meshes representing the different geometries, unstructured volumetric meshes have been created using an octree approach with ICEM CFD 13.0 (ANSYS, Canonsburg, PA, USA). The considered geometries have a typical radius of $2~\text{mm}$ , and the resulting meshes are composed of tetrahedral elements with a side length of 0.24 mm and, at the wall, three layers of prism elements with a total height of $0.08~\text{mm}$ and a side length of $0.12~\text{mm}$ . This mesh configuration was selected according to the mesh convergence study performed in the project @neurist (Aneurist Consortium 2010; Villa-Uriol et al. Reference Villa-Uriol, Berti, Hose, Marzo, Chiarini, Penrose, Pozo, Schmidt, Singh and Lycett2011).

The flow fields have been simulated with the vertex-centred finite volume solver CFX 13.0 (ANSYS) for Navier–Stokes equations, using a second-order advection scheme, a second-order backward Euler transient scheme for unsteady simulations and CFX’ automatic time scale control for steady-state simulations. An incompressible Newtonian fluid with viscosity $\unicode[STIX]{x1D707}=3.5~\text{mPa s}$ and density $\unicode[STIX]{x1D70C}=1066~\text{kg}~\text{m}^{-3}$ (typical values for blood) was considered.

Rigid walls, volumetric inflow rate and parabolic velocity profile in the inlet, with moderate Reynolds numbers ( $\mathit{Re}\leqslant 500$ ), have been imposed as boundary conditions. Steady flows have been simulated for the straight (i) and the double-bend tubes (iii), with a maximum residual of $10^{-5}$ .

Unsteady cyclic flow has been simulated for the arteries with aneurysm (iv). The specific boundary conditions at the inlet (flow rate) and outlets (pressure) were computed for a cardiac cycle using a one-dimensional model of the human vascular tree (Reymond et al. Reference Reymond, Merenda, Perren, Rüfenacht and Stergiopulos2009). The cardiac cycle, with a period of 0.8 s, has been discretized in 160 uniformly distributed time steps ( $5~\text{ms}$ ) and, to reduce the effect of initial transients, the third of three simulated cardiac cycles has been analysed. A maximum residual of $5\times 10^{-3}$ was set as convergence criterion. This corresponds to the state-of-the art protocol for haemodynamic simulation studies defined in the @neurist project. Although not strictly converged in some cases, it represents a balance between accuracy and computational time, being sufficient to illustrate how the ICI could be used in practical applications. The protocol and convergence criterion is stricter for the idealized flow simulations.

For cases (ii) and (v) the flow fields are obtained analytically.

4.3 Streamlines computation

From the flow velocity fields, the streamlines have been obtained by numerical integration using a fourth-order Runge–Kutta algorithm, implemented in the Visualization ToolKit library (Schroeder, Martin & Lorensen Reference Schroeder, Martin and Lorensen2003; VTK community 2014). The integration is controlled by the length step, which we have set to of $1/20$ of the cell length. Inside each cell the velocity is linearly interpolated from the values at the cell nodes. The integration is terminated if either the number of steps reach $10^{5}$ or the speed is below $10^{-100}~\text{m}~\text{s}^{-1}$ . Observe that this termination speed will be virtually never reached, since those stagnant streamlines have null probability according to the natural distribution of streamlines.

The seed points are selected according to the natural distribution (2.8) in either inlet or outlet. The streamlines are propagated in both directions, forwards and backwards. In the inlet case, recirculation is controlled by discarding any seed point for which the backwards propagated streamline crosses again the inlet. In the outlet case, the seed is discarded if the forwards propagated streamline crosses again the outlet. This algorithm implements definition 3 of inlet and outlet.

The streamlines sample is generated on the fly from the flow field while computing the MI. In order to reduce the memory consumption, they are generated in groups of 2k. Then, the inlet and outlet crossings, $\boldsymbol{x}_{In}(\unicode[STIX]{x1D6FE})$ and $\boldsymbol{x}_{Out}(\unicode[STIX]{x1D6FE})$ , are identified and the streamlines are removed to free memory. For the different experiments, MI has been computed from sets of streamlines going from 5k to 100k.

4.4 ICI estimation

The joint probability density (3.2) can be estimated from a sample of $N$ streamlines, $\{\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D6FC}}\}$ , generated from any flow cross-section (in particular, either inlet or outlet) with the natural probability density (2.8):

The factor $p(\unicode[STIX]{x1D6FE})$ disappears because it is taken into account when generating the sample of streamlines. Then, Monte Carlo integration (Robert & Casella Reference Robert and Casella2013) is applied for the computation of the MI in (3.3), using uniform samplings of $M$ points in both the inlet and the outlet. The same procedure is applied for the estimation of $I_{s}(\boldsymbol{X},\boldsymbol{Y})$ , $I_{s}(\boldsymbol{X},\boldsymbol{X})$ and $I_{s}(\boldsymbol{Y},\boldsymbol{Y})$ . From them, $\mathit{ICI}_{s}$ is computed as defined in (3.8). For computation and memory efficiency, the PSF is truncated for distances ${>}3.717\unicode[STIX]{x1D70E}$ , corresponding to the 99.9 % confidence region, and sparse matrices are used in the algorithm.

This algorithm is linear in the number of streamlines, $N$ , and quadratic in the number of points, $M$ . The complete algorithm, including the streamlines computation, has been implemented in C++. All the experiments have been run in a MacBook Pro, 2.8 GHz dual-core Intel Core i7, 16 GB SDRAM. The mean computational time for the full process, from flow field to $\mathit{ICI}_{s}$ , for the scale sequence $s=1/3,1/4,\ldots ,1/20$ , with $N=20\text{k}$ and $M=5\text{k}$ , was 5 min. But this time depends on the geometry extent and the flow complexity.

5.1 Algorithm accuracy and precision evaluation

To check the correctness and accuracy of the algorithm described above for the estimation of the ICI, we have considered the steady flow in a cylindrical straight tube (i) and the pulsatile flow in a subset of the dataset of anatomically realistic vasculatures with aneurysm (iv).

For case (i), the MI can be semi-analytically computed. This computation and the resulting values for $I_{s}(\boldsymbol{X},\boldsymbol{Y})$ are presented in appendix A. Alternatively, the corresponding MI has been also numerically computed from the flow simulated for a tube of radius $R=2~\text{mm}$ and length $L=0.3~\text{m}$ , and a flow rate $Q=5~\text{ml}~\text{s}^{-1}$ , corresponding to $\mathit{Re}=485$ . $I_{s}(\boldsymbol{X},\boldsymbol{Y})$ has been estimated for a series of streamline samples of sizes $N=5\text{k},10\text{k},20\text{k},50\text{k}$ , and of point samples of sizes $M=2\text{k},5\text{k},10\text{k},20\text{k}$ . Each combination has been instantiated 8 times to evaluate also their variability. The resulting values have been compared with the semi-analytical computation, considered as the ground truth (GT). Figure 4 presents the distribution of signed relative errors, $(I_{s}-I_{s}^{(GT)})/I_{s}^{(GT)}$ , obtained with the different sample sizes, for scales $s=1/3,1/10,1/20$ .

Figure 4.

Evaluation of the accuracy and precision of the algorithm implementation for the computation of the MI. For a stationary flow in a cylindrical straight tube, the numerical estimation is compared with the semi-analytical computation. The signed error is evaluated for 3 scales ( $s$ ), 4 streamline sample sizes ( $N$ ) and 4 point sample sizes ( $M$ ).

To estimate the algorithm precision and accuracy in more complex geometries and flows, we have selected a random subset of 10 cases from the dataset of cerebral vasculatures with aneurysms (iv), considering their $\mathit{ICI}_{s}$ computed from the streamlines at one cardiac phase (peak systole). Since no analytical ground truth is available in this case, each $\mathit{ICI}_{s}$ has been computed using $N=100\text{k}$ streamlines and $M=50\text{k}$ points, which is in the limit of the feasible computation, needing for each case 20 h of computation with peaks of 15 GB of memory usage. The obtained values have been taken as ground truth. The ICI has been also computed for $N=5\text{k},10\text{k},20\text{k},50\text{k}$ and $M=2\text{k},5\text{k},10\text{k},20\text{k}$ , instantiated twice for each combination. For each scale, the signed error is given by the difference with respect to the ground truth, $\unicode[STIX]{x1D700}_{s}=\mathit{ICI}_{s}-\mathit{ICI}_{s}^{(GT)}$ . Since the ICI is normalized between 0 and 1, this absolute error is directly interpretable as relative to the index range and comparable between scales. Figure 5 presents the obtained distribution of errors. Table 2 presents the estimated accuracy and precision for 3 different settings for $N$ and $M$ , where the accuracy is quantified by the mean signed error, $\overline{\unicode[STIX]{x1D700}}$ , and the precision by twice its standard deviation, $2\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D700}}$ , representing approximately the 95 % confidence interval ( $\overline{\unicode[STIX]{x1D700}}\pm 2\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D700}}$ ).

Figure 5.

Evaluation of the accuracy and precision of $\mathit{ICI}_{s}$ computation on the flow at peak systole in 10 vasculatures with aneurysm. For each scale ( $s$ ), the signed error against the ground truth is evaluated for several streamline sample sizes ( $N$ ) and point sample sizes ( $M$ ).

Table 2.

Accuracy ( $\overline{\unicode[STIX]{x1D700}}$ ) and precision ( $2\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D700}}$ ) of $\mathit{ICI}_{s}$ for 3 scales and 3 settings of $N$ and $M$ , estimated from the flow in 10 vasculatures with aneurysm. The values are expressed as $\overline{\unicode[STIX]{x1D700}}\pm 2\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D700}}$ (approximate 95 % confidence interval).

As expected, both accuracy and precision increase with the number of streamlines ( $N$ ) and with the number of points ( $M$ ). For the straight tube, the convergence to the semi-analytical values evidences the correctness of the algorithm implementation.

The accuracy is degraded from coarser scales (large $s$ ) to finer scales (small $s$ ), and from larger samples sizes, $N$ and $M$ , to smaller ones. This bias for small samples is a known effect for experimental MI estimators (Paninski Reference Paninski2003). It is interesting that the bias introduced by small $N$ and $M$ acts in opposite directions. Reducing the number of streamlines ( $N$ ) should presumably reduce the apparent complexity, which explains the observed decrease in the ICI (figure 5), specially for $s=1/20$ . Accordingly, this bias dependence with $N$ is not observed for the simple flow in the straight tube (figure 4). For smaller $s$ , the area covered by the PSF is smaller, requiring denser samples. This explains the bias increase when $s$ decreases, for the same $N$ and $M$ .

The estimation of ICI is most precise for the intermediate scale, reducing the precision towards both extremes. Two competing effects appear here. The variability of the MI estimator increases for small sample sizes compared with the observation scale. Thus, for the same $N$ and $M$ , the MI variability decreases when $s$ increases. However, the ICI normalization factor also decreases, magnifying the ICI errors for the same MI errors.

From the obtained values, a reasonable selection seems to be $N=20\text{k}$ and $M=5\text{k}$ , involving an acceptable mean computational time of 5 min per case, representing a balance between accuracy and computational time and memory. Hence, these algorithm settings have been used in the subsequent experiments.

5.2 Streamline distribution invariance

Although the natural probability distribution of streamlines is analytically invariant to the cross-section generating them, this invariance can be broken by the numerical approximations involved in the flow simulations. As shown in § 2.3, the source of variation of the streamline distribution is the integral of the flow divergence, $\unicode[STIX]{x1D709}(\unicode[STIX]{x1D6FE})$ , along each streamline, $\unicode[STIX]{x1D6FE}$ , as defined in (2.18).

For testing this invariance and the influence of the numerical errors, we have considered the unsteady flow cases (iv) in cerebral vasculatures. For each case, 5 different streamline samples have been generated either forward from the inlet or backward from the outlet, and with either the natural or homogeneous distribution.

(FNa)

Forward, natural distribution (Sample A).

(FNb)

Forward, natural distribution (Sample B).

(BN)

Backward, natural distribution.

(FH)

Forward, homogeneous distribution.

(BH)

Backward, homogeneous distribution.

Here, homogeneous denotes the distribution of streamlines obtained from a homogeneous sample of seed points at the inlet (forward) or outlet (backwards). It is used as a naive alternative against which to compare the natural distribution. The corresponding MI has been computed from each of these samples. The flow divergence has been also computed and integrated along each streamline, providing $\unicode[STIX]{x1D709}(\unicode[STIX]{x1D6FE})$ .

Figure 6 presents the MI difference, $\unicode[STIX]{x0394}I$ , between two streamlines samples (inter-sample variability) for 3 different pairs of samples: (FNa)–(FNb), (FNa)–(BN) and (FH)–(BH). The MI difference is plotted against the third quartile of the absolute value of the divergence integrals, $\tilde{\unicode[STIX]{x1D709}}=\text{Q}_{3}(|\unicode[STIX]{x1D709}(\unicode[STIX]{x1D6FE})|)$ . We have considered the linear model $\unicode[STIX]{x0394}I=a+b\tilde{\unicode[STIX]{x1D709}}+\unicode[STIX]{x1D700}+\unicode[STIX]{x1D6FD}\tilde{\unicode[STIX]{x1D709}}$ , with two independent normally distributed errors: a constant noise, $\unicode[STIX]{x1D700}\sim {\mathcal{N}}(0,\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D700}})$ , and an effect linearly increasing with the divergence integral, $\unicode[STIX]{x1D6FD}\tilde{\unicode[STIX]{x1D709}}\sim {\mathcal{N}}(0,\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FD}}\tilde{\unicode[STIX]{x1D709}})$ . The fitted line and 95 % confidence intervals for the differences, are given by the minimum likelihood estimates of the linear model parameter, $a$ , $b$ , $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D700}}$ , $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FD}}$ . The statistical significance of the dependency of the error with the divergence integral ( $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FD}}>0$ ) has been further confirmed by the studentized Breusch–Pagan test (Koenker Reference Koenker1981).

Figure 6.

Dependence of the inter-sample variability of the MI on the divergence integral. The differences are plotted for two streamline samples generated at the inlet (FNa–FNb), sample at the inlet minus sample at the outlet generated by the natural distribution of streamlines (FNa–BN) and generated by homogeneous distribution (FH–BH).

As expected, the difference, $\unicode[STIX]{x0394}I$ , for two samples equally generated from the inlet, (FNa)–(FNb), shows no bias and the 95 % confidence interval is independent of the divergence integral.

For (FNa)–(BN), $\unicode[STIX]{x0394}I$ shows a bias and a 95 % confidence interval significantly dependent on the divergence integral ( $p=6\times 10^{-6}$ ). But the limit for null divergence coincides with the differences obtained from (FNa)–(FNb). This is in agreement with the expected invariance of the natural distribution of streamlines with respect to the generating cross-section.

In contrast, for (FH)–(BH), the 95 % confidence interval is independent of $\tilde{\unicode[STIX]{x1D709}}$ and significantly larger than the confidence interval for (FNa)–(FNb). This reflects the importance of using the natural distribution of streamlines.

6.1 Rotating concentric cylinders with axial flow

The flow resulting in the annulus between two concentric cylinders, of inner radius $R_{1}$ and outer radius $R_{2}$ , rotating with different angular velocities, $\unicode[STIX]{x1D714}_{1}$ and $\unicode[STIX]{x1D714}_{2}$ respectively, is a well-studied case with a long history (Childs Reference Childs2010). In general, the flow can be unstable involving Taylor vortices. For low angular velocity or high viscosity, with Reynolds and Taylor numbers below a critical value, the steady flow is cylindrically symmetric and can be expressed analytically. For the case with static outer cylinder, $\unicode[STIX]{x1D714}_{2}=0$ , and a constant pressure gradient in the axial direction, $\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}z=-p^{\prime }$ , a steady incompressible Newtonian flow under no-slip condition can be expressed in cylindrical coordinates as

where $k=R_{1}/R_{2}$ and $\hat{r}=r/R_{2}$ . This velocity field can be integrated to obtain the streamlines:

Thus, for inlet and outlet defined by the cross-sections at $z=z_{In}$ and $z=z_{In}+L$ , a streamline crossing the inlet at $(r_{In},\unicode[STIX]{x1D703}_{In},z_{In})$ , will cross the outlet at $(r_{Out},\unicode[STIX]{x1D703}_{Out},z_{Out})$ , where

The local flux at the inlet or outlet is given by $\unicode[STIX]{x1D711}_{\boldsymbol{v}}^{{\mathcal{S}}}=v^{z}$ , which gives the total volumetric flow rate

and $A=\unicode[STIX]{x03C0}({R_{2}}^{2}-{R_{1}}^{2})$ is the cross-sectional area. After normalization we obtain the natural probability density of streamlines

Finally, expanding $v^{\unicode[STIX]{x1D703}}$ and $v^{z}$ in (6.3), we get

Observe that $\unicode[STIX]{x1D706}$ is a dimensionless quantity, which can be interpreted as the product of two ratios:

where the Taylor number (or rotational Reynolds number), $\mathit{Ta}$ , and the (axial) Reynolds number, $\mathit{Re}$ , are defined as in Lueptow, Docter & Min (Reference Lueptow, Docter and Min1992) and Wereley & Lueptow (Reference Wereley and Lueptow1999), so that their ratio is the quotient of the velocity of the interior cylinder, $v_{1}^{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D714}_{1}R_{1}$ , by the bulk axial velocity, $\overline{v}^{z}=Q/A$ .

We have run experiments for three values of the ratio between internal and external radius, $k=0.1$ , 0.5 and 0.8, and an exponential sequence of values for $\unicode[STIX]{x1D706}$ , ranging along several orders of magnitude from 0.01 to 100. For each value of $k$ and $\unicode[STIX]{x1D706}$ , a sample of crossing points at the inlet and outlet from $20$ k streamlines have been generated following (6.5) and (6.6). The corresponding ICI has then been computed from each of these samples.

Figure 7 shows the evolution of ICI when $\unicode[STIX]{x1D706}$ increases, for two scales, $s=1/3$ and $1/20$ , and for the three radius ratios. Different combinations of fluid properties or boundary conditions can result in the same value of $\unicode[STIX]{x1D706}$ . Two simple interpretations of the curves are that they represent either the ICI evolution when increasing the distance, $L$ , from inlet to outlet, or the ICI change when increasing the rotational velocity, $\unicode[STIX]{x1D714}_{1}$ .

Figure 7 also illustrates the degree of mixing by rendering how the azimuthally coloured inlet is propagated by the streamlines, producing deformed images by their intersection at the outlet. This gives an intuition of the mixing corresponding to each value of the ICI.

The ICI increases faster and further for coarser scales. In particular, $\mathit{ICI}_{1/3}$ converges very close to 1 for large $\unicode[STIX]{x1D706}$ , while $\mathit{ICI}_{1/20}$ is always below 0.9. At first glance, we could expect that $\mathit{ICI}_{s}$ converges to 1 for any scale. However, this laminar rotating flow produces no mix in the radial direction. Thus, there is always some amount of information preserved. The difference between both scales is larger for the smallest radius ratio, $k=0.1$ , while for $k=0.8$ the ICI of both scales approach each other. The value of $\mathit{ICI}_{s}$ also increases with $k$ , since then the contribution of the radial direction is reduced with respect to the angular one.

Figure 7.

ICI and illustration of the mixing produced by a steady flow in annulus with rotating inner cylinder and axial pressure gradient, for 3 radius ratios. The ICI at two scales ( $s$ ) is plotted as function of the dimensionless parameter $\unicode[STIX]{x1D706}$ . For each radius ratio, the azimuthally coloured inlet (left column) is propagated by the streamlines, producing deformed images by their intersection at the outlet for $\unicode[STIX]{x1D706}=0.01$ , $0.1$ , $1$ , $10$ and $100$ .

6.2 Double-bend tubes

One of the simplest 3-D geometries is a tube including the sequence of two non-coplanar $90^{\circ }$ bends. Lee et al. (Reference Lee, Parker, Caro and Sherwin2008) investigated the characteristics of the steady flow obtained under 5 different inter-bend angles ( $\unicode[STIX]{x1D6FC}=0^{\circ },45^{\circ },90^{\circ },135^{\circ },180^{\circ }$ ) and 2 Reynolds numbers ( $\mathit{Re}=125,500$ ), with tube radius $R=2~\text{mm}$ and bend radius $R_{bend}=8~\text{mm}$ . We have reproduced the same geometries (figure 8) and flow conditions considered there, investigating the ICI obtained in each case.

For each angle, a series of cross-sections have been generated every $1~\text{mm}$ along the tube centreline, numbered according to their position relative to the inter-bend centre (see figure 8). The ICI has been computed from the inlet cross-section at $-30~\text{mm}$ , to each of the outlets in the cross-section series, from position $-20~\text{mm}$ to $40~\text{mm}$ . Figure 9 shows the ICI obtained along this sequence of outlets, for the 5 angles and 2 Reynolds numbers and for 2 scales: $s=1/3$ and $1/20$ .

Figure 8.

Double-bend tubes with inter-bend angles $0^{\circ }$ , $45^{\circ }$ , $90^{\circ }$ , $135^{\circ }$ and $180^{\circ }$ . For each angle, a series of cross-sections have been generated along the tube centreline, each separated 1 mm from the previous one. They are numbered according to their position relative to the inter-bend centre. The inter-bend centre ( $0$ ), the first ring at the beginning of the first bend ( $-12$ ) and the last ring at the end of the second bend ( $12$ ) are shown in yellow.

The inlet has been coloured radially from red in the centre to blue in the perimeter, adding 6 white radius. This image was propagated by the flow streamlines, producing deformed images by their intersection with each cross-section. These images illustrates the degree of mixing and complexity of the flow from the inlet to each of the series of outlets. They are presented in figure 9 at positions $-$ 5, 0, 12 and 40 mm, for each angle and Reynolds, together with the corresponding ICI.

Figure 9.

ICI and illustration of the mixing along double-bend tubes with 5 different inter-bend angles, and for two Reynolds numbers. The corresponding $\mathit{ICI}_{s}$ is computed for two scales, $s=1/3$ (circles) and $s=1/20$ (triangles), from the inlet to each of the outlets in the cross-section series (see figure 8). To visualize the degree of mixing, the radially coloured inlet is propagated by the streamlines, producing deformed images by their intersection with each cross-section. They are presented at positions $-5$ , $0$ , $12$ and $40$  mm, for each angle and Reynolds number.

When comparing the 5 angles, both the coloured cross-sections and the ICI coincide for outlets along the first bend and the initial portion of the second bend, where the dynamics is essentially characterized by the Dean vortices developed by the first bend. The divergence between angles starts earlier for $\mathit{Re}=500$ , as seen in both scales. This can be explained since the Dean vortices generated by the second bend are stronger and faster developed for higher Reynolds number.

The ICI for the planar bend with $\unicode[STIX]{x1D6FC}=180^{\circ }$ is almost the same for both Reynolds numbers at the end of the second bend ( $12~\text{mm}$ ). This can be understood since the higher energy in the Dean vortices for $\mathit{Re}=500$ is compensated by the higher axial velocity, so that the streamlines are similar. Lee et al. (Reference Lee, Parker, Caro and Sherwin2008) mention that at $\mathit{Re}=500$ the strength of the primary Dean vortices is sufficient to lead to secondary and even tertiary vortices. They do not explicitly state this not happening at $\mathit{Re}=125$ , but any possible differences between both does not result in a clear difference in ICI. However, the ICI still grows significantly after the second bend for $\mathit{Re}=500$ , indicating that the vortices are more persistent. This effect is common also for all the angles, so that the ICI at $40~\text{mm}$ is always larger for $\mathit{Re}=500$ , and not even completely developed, in contrast with $\mathit{Re}=125$ .

As expected, the ICI at the outlet is higher for non-planar than for planar models. Both $\unicode[STIX]{x1D6FC}=45^{\circ }$ and $\unicode[STIX]{x1D6FC}=90^{\circ }$ present very similar coloured cross-sections and ICI. They share the highest ICI in both scales and both Reynolds numbers. This is in agreement with the observation in Lee et al. (Reference Lee, Parker, Caro and Sherwin2008) that there appears to be an optimal inter-bend angle between $45^{\circ }$ and $90^{\circ }$ to achieve asymmetry in the Dean vortex pattern.

In contrast with $\unicode[STIX]{x1D6FC}=180^{\circ }$ , for $\unicode[STIX]{x1D6FC}=0^{\circ }$ the ICI is substantially higher at $\mathit{Re}=500$ than at $\mathit{Re}=125$ . Apparently the reverse Dean vortices developed in the second bend, partially undo the mixing for low Reynolds number. This even produces a reduction of ICI for $\mathit{Re}=125$ at the end of the second bend. This effect can be seen in the coloured cross-sections, with the reduction of the deformation in the centre at $12~\text{mm}$ .

For $\mathit{Re}=125$ , we can also observe a small reduction of the ICI at the end of the first bend, only for scale $s=1/3$ , but shared by all the angles. We have not found a satisfactory explanation of this effect.

7.1 Cerebral arteries with aneurysms

As a case study on non-idealized geometries, we have investigated the ICI obtained for unsteady flows in anatomically realistic cerebral arteries with aneurysms. Aneurysms are pathologic dilations of the vessel wall, and their flow characteristics have been extensively investigated as mechanically and biologically related with the initiation, growth and eventual rupture of aneurysms (Cebral et al. Reference Cebral, Mut, Weir and Putman2011; Geers et al. Reference Geers, Larrabide, Radaelli, Bogunovic, Kim, van Andel, Majoie, VanBavel and Frangi2011; Gopalakrishnan, Pier & Biesheuvel Reference Gopalakrishnan, Pier and Biesheuvel2014; Schnell et al. Reference Schnell, Ansari, Vakil, Wasielewski, Carr, Hurley, Bendok, Batjer, Carroll and Carr2014; Arzani et al. Reference Arzani, Gambaruto, Chen and Shadden2016).

A set of 49 geometries obtained from segmentation of 3-D rotational angiography (3DRA) images including an aneurysm in the middle cerebral artery have been selected from the @neurIST database (Aneurist Consortium 2010; Villa-Uriol et al. Reference Villa-Uriol, Berti, Hose, Marzo, Chiarini, Penrose, Pozo, Schmidt, Singh and Lycett2011). Surface models of patient-specific vascular geometries were obtained with the geodesic active region segmentation method (Bogunović et al. Reference Bogunović, Pozo, Villa-Uriol, Majoie, van den Berg, Gratama van Andel, Macho, Blasco, San Román and Frangi2011) and manually corrected for topological errors using the software suite @neuFuse (Villa-Uriol et al. Reference Villa-Uriol, Berti, Hose, Marzo, Chiarini, Penrose, Pozo, Schmidt, Singh and Lycett2011). The vasculature of interest includes at least 12 vessel diameters upstream and 4 vessel diameters downstream from the aneurysm. Inlets and outlets have been clipped with planes perpendicular to the vessel centreline. Unsteady flow simulations have been performed imposing flow rate (inlet) and pressure waveforms (outlets) extracted from a one-dimensional model of all the large arteries in the human body (Reymond et al. Reference Reymond, Merenda, Perren, Rüfenacht and Stergiopulos2009), for a cardiac cycle of period $0.8~\text{s}$ .

For the resulting unsteady incompressible flow, both versions of ICI could be applied: considering either streamlines or world-lines. We have considered here the streamlines for two physiologically relevant cardiac phases: peak systole (PS) and end diastole (ED). This allows us to directly compare the results to some existing subjective classifications of aneurysm flow complexity and stability, defined also for these two cardiac phases.

The considered control volume, ${\mathcal{D}}_{{\mathcal{S}}}$ , comprises the aneurysm and a portion of each adjacent vessel. We have selected inlet and outlet at one vessel diameter from the aneurysm neck, as illustrated in figure 10. This is the same criterion introduced by Millan et al. (Reference Millan, Dempere-Marco, Pozo, Cebral and Frangi2007) for aneurysm morphological characterization. The limiting cross-sections were automatically computed based on the aneurysm removal methodology proposed by Ford et al. (Reference Ford, Hoi, Piccinelli, Antiga and Steinman2009) and implemented in the Vascular Modelling Toolkit (Piccinelli et al. Reference Piccinelli, Veneziani, Steinman, Remuzzi and Antiga2009). Using vessel centrelines and Voronoi diagrams, the aneurysm neck and vessel bifurcations were detected and used as a reference to reliably identify the cross-sections.

For this control volume, the outlet consists of multiple, separate regions. But this does not affect the definition of ICI.

Figure 10.

Example of the cross-sections automatically selected for the inlet vessel and each of the two outlet vessels, at one diameter from the aneurysm neck at the bifurcation.

7.2 ICI as a function of scale

For each case, $\mathit{ICI}_{s}$ was obtained for the sequence of scales, $s=1/3,1/4,\ldots ,1/20$ . This provides ICI as a function of the scale. Figure 11, displays this function for several representative examples.

Figure 11.

ICI as function of scale and illustration of mixing in typical cases from the aneurysm dataset. The line-connected symbols are the ICI at PS of examples representing the minimum, maximum, median and first and third quartiles, for $s=1/20$ . The case with largest change in the order from $s=1/3$ to $s=1/20$ pass from the first to the fourth quartile. The ICI mean, and the variation explained by the first 2 principal components are also plotted. Each of the outlets have been coloured by two different patterns: by homogenous colour identifying each outlet vessel as blue or red, and radially with a rainbow colour map. Their propagation backwards by the streamlines produces respectively two coloured cross-sections at the inlet.

In general, $\mathit{ICI}_{s}$ increases with the scale $s$ . The obtained continuous and smooth behaviour makes evident a high correlation between the values at different scales. Thus, a dimensionality reduction can be obtained, for instance, by principal component analysis. Effectively, the first principal component (PC) explains 99.1 % of the variance, and the 2 first PCs explain 99.9 % of the variance. These two PCs are also plotted in figure 11.

An illustration of the mixing corresponding to each of the ICI values is also presented by coloured cross-sections, similar to the ones in previous cases, but adapted to the multiplicity of the outlet. Each of the outlets have been coloured by two different patterns: radially from red to blue, and by homogenous colour identifying each outlet vessel. Two coloured cross-section are respectively obtained at the inlet, by their propagation backwards along the streamlines.

7.3 Comparison with subjective assessment of flow complexity

For each aneurysm, the flow field has been qualitatively assessed by an expert, through visual inspection of the streamlines, according to four qualitative variables (Cebral et al. Reference Cebral, Mut, Weir and Putman2011; Geers et al. Reference Geers, Larrabide, Radaelli, Bogunovic, Kim, van Andel, Majoie, VanBavel and Frangi2011). Three variables are measured at PS: flow complexity (simple or complex), inflow concentration (concentrated or diffuse) and impingement size (small or large). The fourth variable, flow stability (stable or unstable), is measured by comparing the flow patterns at PS and ED.

We have investigated the correlation of the proposed quantitative ICI at PS with the subjective qualitative classification. For each qualitative variable, the corresponding two categories have been plotted as a box plot of $\mathit{ICI}_{1/10}$ (figure 12). Non-parametric Mann–Whitney $U$ test (Mann & Whitney Reference Mann and Whitney1947; Hollander, Wolfe & Chicken Reference Hollander, Wolfe and Chicken2013) has been applied to each pair of categories. Three of them gave statistically significant differences (flow complexity: $p=9\times 10^{-3}$ , flow stability: $p=2\times 10^{-4}$ , impingement size: $p=4\times 10^{-4}$ ). Only inflow concentration did not show significant differences. Surprisingly, flow stability showed significant differences for the ICI exclusively measured at PS, although it is supposed to account for the variation between ED and PS. Probably, this reflects the fact that more complex flows are also more unstable.

Figure 12.

Correlation of ICI with each of the 4 subjective flow classifications. The significant difference between each pair of categories is assessed by non-parametric Mann–Whitney $U$ test. $^{\ast \ast }(0.01\geqslant p>0.001)$ , $^{\ast \ast \ast }(0.001\geqslant p>0.0001)$ .

7.4 Comparison between peak systole and end diastole

Figure 13 presents a dot plot with $\mathit{ICI}_{1/10}$ obtained for each aneurysm at PS and ED. The values spread quite homogeneously across the ICI range. In general, the ICI obtained for PS is larger than the one obtained for ED, but the opposite behaviour is also observed.

Figure 13.

Comparison between ICI of the flow at peak systole and end diastole for each aneurysm, at scale $s=1/10$ . The ICI is generally larger at PS than at ED, but the opposite behaviour is also observed. The aneurysms have been ordered according to their ICI at ED.

The ICI variation from ED to PS is presented in figure 14, separated according to the qualitative classification into stable and unstable. A clear correlation appears, which is confirmed by one-tailed Mann–Whitney $U$ test, giving a high significant difference ( $p=5\times 10^{-4}$ ). For completeness, Mann–Whitney test was also applied to the ICI variation for the classes obtained from the qualitative flow complexity, inflow concentration and impingement size. None of them provided significant differences.

Figure 14.

Subjective flow stability classes against absolute value of the difference of ICI between peak systole and end diastole. A strong significant difference ( $p=0.0005$ , one-tailed Mann–Whitney $U$ test) is found between subjectively stable and unstable flows.

8.1 Two-dimensional periodic cellular flows

Periodic cellular flows are simply defined models that can, however, show very complex evolutions. We have considered a modification of the 2-D unsteady model used by Zu, Chen & Xin (Reference Zu, Chen and Xin2015):

defined in Cartesian coordinates $\boldsymbol{x}=(x,y)$ for the square domain $\unicode[STIX]{x1D6FA}=[0,1]\times [0,1]$ with periodic boundary conditions (torus). The base flow ( $\unicode[STIX]{x1D70C}=1,\unicode[STIX]{x1D703}=0$ ) produce two counter-rotating cells centred at $(1/4,3/4)$ and $(3/4,1/4)$ , and two saddles at $(1/4,1/4)$ and $(3/4,3/4)$ . The parameter $\unicode[STIX]{x1D703}$ introduce a time-dependent perturbation, making it an unsteady flow. The original model did not include the density $\unicode[STIX]{x1D70C}$ . We have considered a density

controlled by the parameter $\unicode[STIX]{x1D705}$ , which introduces an inhomogeneity, making the flow compressible. This density is maximum at the saddles and minimum at the rotating centres. The resulting unsteady compressible flow satisfies the continuity equation (2.28). Thus, the natural distribution of world-lines is invariant, and we can consider the corresponding version of ICI.

For this flow, we consider the spatio-temporal control volume ${\mathcal{D}}_{{\mathcal{S}}}=\unicode[STIX]{x1D6FA}\times {\mathcal{T}}$ given by the full spatial cyclic domain $\unicode[STIX]{x1D6FA}$ (torus) along a lapse of time ${\mathcal{T}}=[0,T]$ . Thus, the inlet and outlet are not spatial surfaces, but the full space, $\unicode[STIX]{x1D6FA}$ , at instances $t=0$ and $t=T$ , respectively.

Figure 15 illustrates the resulting flow for different values of $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D705}$ , and for different time lapses. The ICI for each of these parameter values and time lapses is displayed for a series of scales from $s=0.01$ to $s=0.5$ . The flow propagation of a vertical colour pattern is also shown for each of the parameter values.

Figure 15.

ICI as a function of scale for the 2-D periodic cellular flow with $\unicode[STIX]{x1D705}=0$ (incompressible) and $\unicode[STIX]{x1D705}=0.9$ (compressible), with 4 different values of $\unicode[STIX]{x1D703}$ (unsteadiness) and for 3 time lapses. The legend is presented at the bottom, together with the illustration of the propagation by each flow of a vertical colour pattern.

The dominant factor is the increase of ICI with the time lapse. The increase of ICI with the parameter $\unicode[STIX]{x1D703}$ is pronounced for all scales and time lapses, reflecting that the higher the unsteadiness, the more complex the flow. The ICI for compressible flows ( $\unicode[STIX]{x1D705}=0.9$ ) is also higher than for incompressible flows ( $\unicode[STIX]{x1D705}=0$ ), although the difference is less marked. It is more evident for the steady flow ( $\unicode[STIX]{x1D703}=0$ ) and intermediate time ( $t=1$ ), for which the difference between steady and unsteady flow is reduced. This is in agreement with the intuitive impression from the propagated colour pattern.

The curves of $\mathit{ICI}_{s}$ increase with the scale. But each of the cases presents a clear maximum slope at different scales, where the complexity increases faster. The plots evidence that the complexity is progressively propagated to smaller scales with time, in agreement with the colour patterns. We can also see that the higher the unsteadiness ( $\unicode[STIX]{x1D703}$ ), the smaller the scale showing fast ICI increase.

8.2 Test of the limit of ICI for infinitesimal scales

In appendix B, we demonstrate (B 13) that $\lim _{s\rightarrow 0}\mathit{ICI}_{s}=0$ , and that the dominant term for infinitesimal scales is proportional to $-1/\log (s)$ .

For the proposed cellular flow, the probability densities at inlet at outlet coincide and are given by the density,

Thus, $\unicode[STIX]{x0394}=0$ , and (B 13) becomes

Observe that the term $\log (|\unicode[STIX]{x1D64C}_{\boldsymbol{x}}|/4)$ shares some elements with the finite-time Lyapunov exponents (Tang & Boozer Reference Tang and Boozer1996), since it can be expressed in terms of the eigenvalues of the Cauchy–Green deformation tensor $\unicode[STIX]{x1D645}_{\boldsymbol{x}}^{\text{T}}\unicode[STIX]{x1D645}_{\boldsymbol{x}}$ . Actually, the higher the Lyapunov exponents, the higher will be $\unicode[STIX]{x1D6F9}$ . However, Lyapunov exponents are intensive quantities (average flow characteristics per unit of time), whereas ICI is extensive (total flow characteristics in a finite spatio-temporal domain). Thus, in contrast with Lyapunov exponents, $\unicode[STIX]{x1D6F9}$ is not normalized by the time lapse, so that it is expected to unboundedly grow with time.

We have numerically computed the factor $\unicode[STIX]{x1D6F9}$ for flows with different values of $\unicode[STIX]{x1D703}=0,1,2,3$ and $\unicode[STIX]{x1D705}=0,0.5,0.9$ , and for different time lapses, $t=0.1,1,10$ . To compute the integral, the factor $\log (|\unicode[STIX]{x1D64C}_{\boldsymbol{x}}|/4)$ has been computed in a regular grid of $100\times 100$ points. For this, the Jacobian matrix, $\unicode[STIX]{x1D645}_{\boldsymbol{x}}$ , has been computed from the flow equations applying the method described in (Sandri Reference Sandri1996), using a fourth-order Runge–Kutta algorithm with time step $10^{-4}$ .

Figure 16 compares the behaviour of $\mathit{ICI}_{s}$ for small scales to the deduced behaviour dependent on $\unicode[STIX]{x1D6F9}$ . The left plot displays $\mathit{ICI}_{s}$ against $-1/\log (s)$ together with the expected limiting straight line of slope $\unicode[STIX]{x1D6F9}$ . The right plot complements it by displaying $-\log (s)\times \mathit{ICI}_{s}$ , whose limit for $s\rightarrow 0$ is directly $\unicode[STIX]{x1D6F9}$ , being more easily observable. In all cases, the expected limit seems to be satisfied. This is more evident for $t=0.1$ and $t=1$ , for which a linear behaviour is reached for the shown scales, $s\geqslant 0.01$ . For $t=10$ , this is not still reached, but the limit, $\unicode[STIX]{x1D6F9}$ , seems compatible with a smooth propagation of $-\log (s)\times \mathit{ICI}_{s}$ . Observe that they show ( $t=1$ ) or suggest ( $t=0.1$ and $t=10$ ) a maximum for some scale. The larger $t$ and $\unicode[STIX]{x1D703}$ , the smaller this scale.

Figure 16.

$\mathit{ICI}_{s}$ for small scales in comparison to the expected leading-order term $\unicode[STIX]{x1D6F9}/\log (s^{-1})$ . They are shown for the periodic cellular flow with $\unicode[STIX]{x1D705}=0.5$ , for 4 values of $\unicode[STIX]{x1D703}$ and 3 time lapses $t$ . (a)  $\mathit{ICI}_{s}$ against $1/\log (s^{-1})$ (connected symbols) and limit slope $\unicode[STIX]{x1D6F9}$ (dashed lines). (b)  $\mathit{ICI}_{s}\times \log (s^{-1})$ against $1/\log (s^{-1})$ (connected symbols) and limit $\lim _{s\rightarrow 0}\mathit{ICI}_{s}\times \log (s^{-1})=\unicode[STIX]{x1D6F9}$ (symbols at $s=0$ ).