There is a growing medical interest in the prediction of the flow and forces inside cerebral aneurysms [24, 52], with the ultimate goal of supporting medical procedures and decisions by presenting viable scenarios for intervention. The clinical background of intracranial aneurysms and subarachnoid hemorrhages is well introduced in the literature such as [46, 51]. These days, with the development of high-precision medical imaging techniques, the geometry and structure of blood vessels and possible aneurysms that have formed, can be accurately determined. To date, surgeons and radiologists had to make decisions about possible treatment of an aneurysm based on size, shape and location criteria alone. In this chapter we focus on the role of Computational Fluid Dynamics (CFD) for therapeutic options in the treatment of aneurysms. The tremendous potential of CFD in this respect was already anticipated in . The value of numerical simulations for treating aneurysms will likely increase further with better quantitative understanding of hemodynamics in cerebral blood flow.
Ultimately, we aim to support the medical decision process via computational modeling. In particular, CFD allows to add qualitative and quantitative characteristics of the blood flow inside the aneurysm to this complex decision process. We propose to compute the precise patient-specific pulsatile flow in all spatial and temporal details, using a so-called ‘Immersed Boundary’ (IB) method. This requires a number of steps, from preparing the raw medical imagery to define the complex patient-specific flow domain, to the execution of high-fidelity simulations and their detailed interpretation in terms of flow visualization and the extraction of quantitative measures of relevance to medical practice. We compute the flow inside the aneurysm to predict high and low stress regions, of relevance to the possible growth of an aneurysm. We also visualize vortical structures in the flow indicating the quality of local blood circulation. We show that, as the size of the aneurysm increases, qualitative transitions in the flow behavior can arise, which express themselves as high-frequency variations in the flow and shear stresses. These variations could quantify the level of risk associated with the growing aneurysm. Such computational modeling may lead to a better understanding of the progressive weakening of the vessel wall and its possible rupture after long time.
In this chapter we present a numerical model for the simulation of blood flow inside cerebral aneurysms. A significant amount of work has been done on simulation of flow in the brain and in the cardiovascular system [2, 5, 13, 18, 24, 38, 40]. As a numerical approach, the finite element method is most commonly used to represent geometries of vessels and the flow of blood through them. Often, the data are obtained from rather coarse biomedical imagery. As a result, the highly complex vessel geometry is defined with some uncertainty, and considerable smoothing and interface approximation need to be included to prepare simulations with a body-fitted approach [2, 6, 11]. As an alternative, the IB method was designed primarily for capturing viscous flow in domains of realistic complexity . In particular, we consider a volume penalization method. In this method, fluid is penalized from entering a solid part of a domain of interest by adding a suitable forcing term to the equations governing the fluid flow . This method is also known as ‘fictitious domain’ method  and physically resembles the Darcy penalty method  or the Navier-Stokes/Brinkman equations describing viscous flow on the scale of individual pores in porous domains . Here, we consider in particular the limit in which the porous domain becomes impenetrable and flow in complex solid domains can be represented. This method is discussed as one of the IB methods in the recent review chapter by , and in the sequel will be referred to as ‘volume penalization immersed boundary method’, a label that was also adopted in .
A primary challenge for any CFD method, whether it is a body-fitted method  or an IB method [17, 32, 39], is to capture the flow near solid-fluid interfaces. In this region the highest velocity gradients may occur, leading to correspondingly highest levels of shear stress, but also potentially highest levels of numerical error. In methods employing body-fitted grids, the quality of predictions is directly linked to the degree to which grid-lines can be orthogonal to the solid-fluid interface and to each other. Also, variation in local mesh sizes and shapes of adjacent grid cells is a factor determining numerical error. The generation of a suitable grid is further complicated as the raw data that define the actual aneurysm geometry often require considerable preprocessing steps before any grid can be obtained. These steps include significant smoothing, segmentation and geometric operations eliminating small side vessels that are felt not to be too important for the flow. On the positive side, the main benefit of a body-fitted approach is that discrete variables are situated also at the solid-fluid interface, which makes implementation of no-slip boundary conditions quite straightforward. Hence, in body-fitted approaches the no-slip property can be accurately imposed, but only on a ‘pre-processed’ smoothed and often somewhat altered geometry [6, 11]. These finite element based approaches can be used to predict the patient’s main flow structures of clinical value as suggested by [5, 6, 15, 19].
Capturing flow near complex shaped solid-fluid interfaces is equally challenging in an IB method, as it is in body-fitted approaches. In the IB approach adopted here, the actual geometry of the aneurysm can be extracted directly from the voxel information in the raw medical imagery, without the need for smoothing of the geometry. Grid generation is no issue for IB methods since the geometry of the flow domain is directly immersed in a Cartesian grid. The location of the solid-fluid interface is known only up to the size of a grid cell, and the shape of the interface is approximated using a ‘staircase’ representation, stemming from the fact that any grid cell is labeled either entirely ‘solid’ or entirely ‘fluid’. Refinements in which a fraction between 0 and 1 of a cell can be fluid-filled  are not taken into consideration here. In fact, the medical imagery from which we start has a spatial resolution that is not too high when small-scale details are concerned. This calls for a systematic assessment of the sensitivity of predictions to uncertainties in the flow domain  incorporating also the effects due to adaptations of the domain by smoothing and interface reconstruction as is considered in higher-order methods . Without relaxing the staircase approximation, the problem of capturing near-interface properties can only be addressed by increasing the spatial resolution. This gives an insight into the error-reduction by systematic grid-refinement for flow in complex geometries.
We illustrate the process of predicting flow and forces based on incompressible Newtonian fluid to characterize blood properties. Non-Newtonian corrections can be readily included, however, these typically lead only to modest changes in the predictions and will hence be omitted here. Pulsatile flow forcing is obtained from the direct measurement of the time-dependent mean flow velocity in a vessel during a cardiac cycle. Transcranial Doppler (TCD) sonography is a non-invasive technique, which can be used for this purpose, allowing to measure cerebral blood flow velocity near the actual aneurysm . We consider a full range of physiologically relevant conditions. Understanding flow patterns inside an aneurysm may help to describe long-term effects such as the likelihood of the growth  or even rupture  of the aneurysm, or the accelerated deterioration of the vessel wall due to low shear stress . Regions of high and low shear stress are often visualized as potential markers for aneurysm growth. High shear stress levels were reported near the ‘neck’ of a saccular aneurysm, and may be relevant during the initiating phase . Low wall shear stress has been reported to have a negative effect on endothelial cells and may be important to local remodeling of an arterial wall and to aneurysm growth . A low wall shear stress may facilitate the growing phase and may trigger the rupture of a cerebral aneurysm by causing degenerative changes in the aneurysm wall. The situation is, however, more complex, as illustrated by the phenomenon of spontaneous stabilization of aneurysms after an initial phase of growth . It is still very much an open issue what the precise relation is between shear stress patterns and general circulation on the one hand, and developing medical risks such as aneurysm rupture, on the other hand. In this complex problem, hemodynamic stimuli are but one of many factors.
Cerebral aneurysms are most often located in or near the Circle of Willis  – the central vessel network for the supply of blood to the human brain. Common risk-areas are at ‘T’ and ‘Y’-shaped junctions in the vessels . Treatment of cerebral aneurysms often involves insertion of coils. This coiling procedure represents considerable risk and uncertainty about the long-term stability of coiled aneurysms [45, 47]. Blood vessels and aneurysms are rather complex by their structure and geometrical shapes. The walls of blood vessels contain several layers of different types of biological cells, which provide elasticity to the vessels and play a role in the compliance . The shape of cerebral aneurysms developing in patients can be inferred by using three-dimensional rotational angiography . In this procedure a part of the brain can be scanned, and aneurysms even of a size less than 3 mm can be depicted [3, 48]. This technique allows a reconstruction of three-dimensional arteries and aneurysms and hence an approximate identification of the blood vessels and parts of the soft tissue in the scanned volume. In the IB approach the domain is characterized by a so-called masking function, which takes the value ‘0’ in the fluid (blood) part and ‘1’ in solid (tissue) parts of the domain. The raw angiography data allows for a simple ‘staircase approximation’ of the solid-fluid interface that defines the vessel and aneurysm shape. Individual voxels in the digital data form the smallest unit of localization of the solid-fluid interface. A computational cell is assigned to be ‘solid’ or ‘fluid’ on the basis of the digital imagery. We will adopt the ‘staircase’ geometry representation in this chapter and do not incorporate any additional smoothing of the geometry or sophisticated reconstruction methods.
For a more complete modeling of the dynamics in the vessel system, flow-structure interaction often plays a role . In that case also parameters and models that characterize, e.g., mechanical properties of arterial tissue, influence of brain tissue and the influence of the cerebrospinal fluid are required. The amplitude of the wall motion in intracranial aneurysms was found to be less than 10% of an artery diameter. Despite the rather modest motion of the vessel, long time effects may accumulate. Even modest movement can affect the vessel walls, which might play a role in possible aneurysm rupture as was hypothesized in . For realistic pulsatile flows some movement of the aneurysm walls was observed during a cardiac cycle . In this chapter we take a first step and restrict to developing the IB approach for rigid geometries. This allows to obtain the main flow characteristics inside relatively large cerebral aneurysms for which the wall movement can be neglected .
The organization of this chapter is as follows. In Section 2 we present the computational model, discuss numerical discretization and introduce the IB method for defining complex vessel and aneurysm geometries. We also describe the process of the reconstruction of the geometry from medical imagery. We illustrate steady flow inside a realistic aneurysm geometry in Section 3 and discuss the reliability of numerical predictions. A pulsatile flow and qualitative impression of the flow and forces distribution inside a realistic aneurysm is presented in Section 4. Concluding remarks are in Section 5.
2. Numerical model of blood flow inside human brain
In this section we present the numerical model for simulation of blood flow inside the human brain. We consider blood as an incompressible Newtonian fluid with a constant viscosity. The flow dynamics is governed by the Navier-Stokes equations, which are presented in Subsection 2.1. The numerical method used to solve these equations in 3D is based on skew-symmetric finite-volume discretization and explicit time-stepping using the method as proposed by . This is discussed in Subsection 2.2, in which we also detail the volume penalizing IB method that is adopted to represent the complex vessel structures through which the blood flows. Finally, in Subsection 2.3, we discuss the method used to obtain the precise geometry of realistic blood vessels.
2.1. The Navier-Stokes equations for an incompressible cerebral flow
There are various approaches to model flow of blood in the human brain. A comprehensive overview is given in . In one approach, the blood is approximated as a Newtonian fluid . More refined models, e.g., the Carreau-Yasuda model, include the shear-thinning behavior of blood and allow to capture non-Newtonian rheology [2, 13, 18]. Under physiological flow conditions in sufficiently large arteries non-Newtonian corrections were found to be quite small [6, 13, 18, 38]. The main flow characteristics appeared to be the same as for a Newtonian fluid at somewhat different stress and velocity levels.
We concentrate on the human brain, in particular on arteries of the Circle of Willis. Typical fine-scale structures in the blood are on the order of 10−6m. A length-scale that characterizes the cross section of a typical cerebral vessel inside the Circle of Willis is on the order of 10−3m . This difference in length-scale of three orders of magnitude motivates to approximate blood as an incompressible Newtonian fluid .
The Navier-Stokes equations provide a full representation of Newtonian fluid mechanics, expressing conservation of mass and momentum. The total physical domain Ω, consists of a fluid part Ωf and a solid part Ωs. The interface between the two will be identified as ∂Ω at which no-slip conditions apply. The governing equations are given in dimensional form by:
Here u∗ is the velocity of the fluid, ρ∗ its mass density, p∗ the pressure and f∗ a forcing term that will play a central role in this chapter as it is used to represent the impenetrability of complex shaped solid vessel walls. We denote variables with a physical dimension by an asterisk and render the formulation dimensionless momentarily. By choosing a reference velocity and reference length we can express a reference time scale as Using as reference density we will use a reference pressure as For the forcing term we select a direct volume penalization in which
where ε∗ is a forcing time scale and H is the masking function: H(x) = 1 if x ∈ Ωs and H(x) = 0 if x ∈ Ωf. We set the reference forcing time scale i.e., much smaller than the reference time scale. Here we introduce the dimensionless forcing parameter ε<< 1.
After choosing all reference parameters we obtain the non-dimensional form of the Navier-Stokes equations:
In this chapter we will consider only stationary, i.e., non moving walls. By adding the forcing to the incompressible momentum equations we formally arrive at the Brinkman equation for flow in a porous medium with permeability related to the parameter ε . Note that, with the inclusion of f as in (3) we arrive at a model for the velocity and pressure fields in the entire domain Ω.
The Reynolds number Re is the only parameter which is required to specify the flow conditions. It quantifies the ratio between the magnitude of the (destabilizing) convective transport and the (stabilizing) viscous processes. It is well known that for relatively low Reynolds numbers flow is laminar and steady , which implies a smooth velocity and pressure field. With increasing Reynolds number the flow can develop more detailed vertical structures, e.g., associated with separated flow near abrupt changes in the shape of a vessel. A further increase in Re usually implies that the flow becomes unsteady and the range of vortices becomes much wider . The range of Reynolds numbers arising in the flow in the Circle of Willis, corresponds to laminar, possibly unsteady flow. This will be discussed in more detail momentarily.
2.2. Numerical method for flow simulation using an immersed boundary approach for representing complex geometries
In this subsection we sketch the numerical method used for the simulation of flow through complex shaped domains. First, we describe the direct numerical simulation approach and specify the volume penalization IB method afterwards.
We employ a staggered allocation of the flow variables (u, p) = (u, v, w, p) as basis for our flow solver . In two dimensions this is sketched in Figure 1, where a primary grid cell with the pressure defined in the center and the Cartesian velocity components at the cell surfaces is presented. The locations at which the velocities and the pressure are stored are referred to as the velocity- and the pressure-points, respectively.
The principles of conservation of mass and momentum as expressed in (4) and (5), form the basis for the discrete computational model that is used for the actual simulations. In the Navier-Stokes equations (4) the rate of change of momentum is obtained from the nonlinear convective flux, the linear viscous flux, the gradient of the pressure and the contribution from the forcing term. These contributions to the total flux each have a particular physical character that needs to be represented properly in the discrete formulation. In particular, the convective flux is skew-symmetric, implying that this flux only contributes to the transport of kinetic energy of the solution in physical space; it does not generate nor dissipate this energy. Likewise, the viscous flux contributes only to dissipation of energy, which has to be strictly maintained in a numerical method. We motivate this in some more detail next.
Starting from the original momentum equation without the forcing term
we are interested in the kinetic energy, given by
where u u is the vector inner product. Note, that in (7) we effectively integrate only over Ωf as u = 0 in Ωs. The evolution of the kinetic energy follows from
In order to obtain the integrand in (8) we multiply the Navier-Stokes equation (6) by u. Integrating by parts we can derive the contribution of each of the fluxes in (6). In fact, the convective and pressure terms do not contribute to the evolution of energy, and we find
where ∇u : ∇u = ∂iuj∂iuj in which we sum over repeated indices. This suggests that the energy of any solution decreases in time because of viscous fluxes only.
A more elaborate derivation can be found in , which departs from the expressions
where energy is written in terms of the function inner product (u, u) as defined in (7). This yields the symmetric expression
Since (∇p, u) = −(p,∇ u) and the skew-symmetry ((u ∇)v,w) = −(v, (u ∇)w) we obtain again (9), i.e., no contribution from pressure and the convective flux.
In a discrete setting the Navier-Stokes equations in matrix-vector notation are written as
where uh is the vector containing the discrete velocity solutions (uh, vh, wh), ph is the discrete pressure, Λ is a matrix with the volumes of the grid cells on its diagonal, C and D are the coefficient matrices corresponding to the discretization of the convective ((u ∇)u) and diffusive (−Δu/Re) operators, respectively. The discretization of the pressure gradient is given by −MT, while the coefficient matrix M itself represents the discretization of the divergence operator, integrated over the control volumes .
The discrete approximation for the kinetic energy can be given using the midpoint rule, as
Similar to the continuous case (11) we compute the evolution of the energy in the discrete model as
For a discrete solution we also require the convective conservation of energy, which implies skew-symmetry of the matrix C of the convective operator: C + CT = 0. The two terms related to the numerical pressure gradient can be rewritten as
where the numerical divergence operator M satisfies equation (13). Thus, pressure terms also cancel and do not influence the stability of the spatial discretization. By comparison with the expression for the energy evolution (9), the second term in the right-hand side of (15) should provide a strict decrease of the energy:
This implies that the coefficient matrix D of the diffusion operator is a positive-definite matrix. Thus, discretely we obtain the same properties for the energy decay as in the continuous case.
In this chapter we employ symmetry preserving finite volume discretization and use central differencing of second order accuracy, which maintains explicitly the skew-symmetry in the discrete equations. Since the energy is preserved under the convective operator the skew symmetric discretization allows to obtain a stable solution on any grid. For proper capturing of the solenoidal property (5) of the velocity field we approximate the gradient operator by the transpose of the numerical divergence operator and use a positive definite discretization of the viscous terms, closely following . The contributions of the convective, viscous and pressure gradient fluxes are integrated in time using a generalization of the explicit second order accurate Adams-Bashforth method. Care is taken of accurately representing the skew-symmetry also in the time-integration. Full incorporation would require an implicit time-stepping, which, however, is computationally too demanding. Instead, time-integration starts from a modification of the leapfrog method  with linear inter/extrapolations of the required ‘off-step’ velocities and an implicit treatment of the incompressibility constraint. Optimization for largest stability region of the resulting scheme yields a particular so-called ‘one-leg’ time-integration method, with a mathematical structure that is akin to the well-known Adams-Bashforth scheme. More details can be found in .
For the total computational domain periodic boundary conditions are applied. A special role is played by the forcing term f in the Navier-Stokes equations (4),which represents the volume penalization accounting for solid objects inside and at the boundaries of the flow domain. The role of the forcing term is to yield an accurate approximation of the no-slip condition at solid boundaries. In conventional computational fluid dynamics such a forcing term is not needed since the flow domain is endowed with a body-fitted grid on which the equations are discretized. The grid-lines in such cases are defined such that they either closely follow the contours of the solid boundaries, or they are (preferably) orthogonal to them. In such a discrete formulation the no-slip boundary condition can be imposed easily. The body-fitted grid is efficient if the fluid domain Ωf is not too complex and does not contain too many separate objects around which the fluid should flow . For considerably more complex flow domains or in case the location of the solid-fluid interface is not perfectly known, as in case of medical imagery, the body-fitted grid approach is limited by the generation of suitable meshes. These should not only align with the solid boundaries, but also be sufficiently smooth near these boundaries to allow an accurate solution in the boundary layers . In our discrete model the forcing term contributes strongly to the stiffness of the equations. When an explicit time-stepping method would be adopted for the forcing term, as is done for the other dynamic contributions, this would result in extremely small time-steps in view of numerical stability. Therefore, the linear forcing term is integrated in time using the implicit Euler scheme .
We use an IB method based on volume penalization  to capture flow in complex aneurysm geometries. A key role in our IB method is played by the so-called ‘masking function’ H, which is a binary function in 3D with values ‘0’ inside the fluid and ‘1’ in the solid parts of the domain. In regions where H = 0 the Navier-Stokes system is solved and thus the fluid domain is defined. In the solid regions H = 1 and the forcing is dominant if the non-dimensional parameter ε is very small. We take ε = 10−10 in the sequel. As a result, the momentum equation (4) reduces to ∂tu ≈ −u/ε if |u| >> ε in the solid domain. Hence, any nonzero u is exponentially sent back to 0 on a time-scale ε. If |u| ≤ ε the forcing is not dominant in the solid, but control over |u| is already obtained, i.e., |u| takes on negligible values in the solid.
A schematic representation in 2D of the masking function for a complex domain is presented in Figure 2, where dark cells correspond to fluid cells with masking function H = 0, while white cells represent the solid part of the computational domain. For realistic 3D cases uncertainties in defining the geometry arise because of the finite spatial resolution of the images of the vessels. We introduced and applied so-called numerical ‘bounding’ solutions in , where next to the basic constructed geometry we consider slightly smaller (‘inner’) and slightly bigger (‘outer’) geometries. These bounding geometries differ from the ‘basic’ geometry by maximally one grid cell in terms of the location of the numerical solid-fluid interface. This implies very tight bounding of the basic geometry especially at high grid resolution. We routinely compute the flow in the basic geometry and in two bounding geometries in order to quantify the associated level of uncertainty. The solution and its main characteristics in the basic geometry appeared to be bounded by the inner and outer solutions.
2.3.Masking. function of realistic cerebral aneurysm, reconstructed from 3DRA data
In this subsection we describe the process of construction of the masking function from medical imagery. The initial set of data was obtained from three-dimensional rotational angiography (3DRA) applied to the brains of patients suffering from a brain aneurysm. The 3DRA scans were provided by the St. Elizabeth Hospital in Tilburg (The Netherlands). We choose one example to illustrate our numerical method.
Before simulation of flow we need to convert the geometry into the masking function. The 3D volume data was first segmented to obtain a vessel structure suitable for flow simulations. After this step small branch vessels were removed and the main vessel with the aneurysm on it was retained . An illustration of the geometry obtained at this stage is presented in Figure 3(a). On the right hand side downstream of the aneurysm we observe the vessel to split. Currently, our approach only allows for a single inflow and a single outflow. Therefore the next steps of the construction are ‘cutting’ a part of the vessel and ‘connecting’ the ‘outflow’ of the vessel to its ‘inflow’ by adding an appropriate, smooth connecting vessel. This way we obtain the desired periodic flow model. Adding an artificial connecting vessel can be avoided by allowing actual inflow and outflow boundary conditions. However, for the flow in the direct vicinity of the aneurysm, we expect the influence of the periodicity assumption relatively ‘far’ from the aneurysm to be negligible and the remaining computational model to be suitable for illustrating the flow in realistic aneurysms. In Figure 3(b) this connecting vessel is represented in black. Several steps in the specification of the masking function can be a source of error in the flow prediction. Special care needs to be taken to limit these errors and to understand the sensitivity of the predictions to details of the connecting vessel, locations of cuts etc. etc. Through systematic numerical simulations these sensitivities can be assessed and the reliability of the predictions quantified.
3. Steady flow simulations in realistic cerebral aneurysm
In this section we present numerical results for the flow inside the realistic aneurysm geometry as shown in Figure 3. We simulate steady flow and present the solution (velocity and pressure) and its gradients (in terms of the shear stress) in Subsection 3.1. Subsequently, we consider the reliability of the predictions by presenting the convergence of pressure drop, velocity and shear stress in Subsection 3.2.
3.1. Motivation and exact definition of the reference case
We simulate flow in the geometry as introduced in Figure 3. First, we need to specify the correct scale and corresponding non-dimensional parameters, to define this application in detail. Subsequently, we will visualize the solution in a number of qualitative ways to appreciate the complexity of the flow structures that develop.
In order to specify the computational model we need to choose reference scales for the length and velocity, and specify the kinematic viscosity. These quantities allow to compute the Reynolds number Re, which is the only parameter that is required for the dimensionless formulation. From literature one may find a range of values characteristic of these reference scales, which implies some degree of freedom in setting the precise values that are physiologically relevant. In particular, the value of the flow velocity is subject to the largest uncertainty when comparing literature, although also the viscosity displays considerable variation. We motivate our choices next to provide a point of reference.
The raw data of the 3DRA scan of the aneurysm consists of a grid of 2563 voxels. The voxel width is 0.1213 mm. This implies that the total physical length of the flow domain is 3.10528 cm. As reference length we take a characteristic scale representative of the average radius of the cerebral vessel in this part of the Circle of Willis. We extract R∗ = 1.94 mm from the 3DRA data. This value is quite similar to  who adopt 2.5 mm, and also consistent with , who suggests a scale of 2.1 ± 0.4 mm. Hence, in the non-dimensional setting we work in a domain of total ‘length’ of 31.0528/1.94 ≈ 16. The computational domain, enclosing the aneurysm geometry is in fact 16×8×16 in the non-dimensional formulation. Simulations will be done at grid resolutions 128 × 64 × 128, 64 × 32 × 64 and 32 × 16 × 32. These resolutions define the grid refinements that we use for the convergence analysis.
Next to the length-scale, also the viscosity and the velocity scales need to be set. From literature we infer that the mass density ρ∗ is in the range 1025 ≤ ρ∗ ≤ 1125 kg/m3, while the dynamic viscosity of blood is reported to be 3 10−6 ≤ μ∗ ≤ 4 10−6 Pa s. Specifically, choosing typical values for the mass density of blood ρ∗ = 1060 kg/m3 and the dynamic viscosity of blood μ∗ = 3.2 10−3 Pa s implies a kinematic viscosity ν∗ = μ∗/ρ∗ = 3.01 10−6 m2/s. Finally, the reference flow-rate as proposed by  and  is 245 ± 65 ml/min, showing an uncertainty of about 25%. This range of values was obtained on the basis of either 3D MR angiograms or TCD measurements. The corresponding range for the velocity scale can be extracted from this as u∗ = Q∗/(π(R∗)2) = 0.345 ± 0.09 m/s. This is consistent with the range 0.34±0.087 m/s as obtained by  on the basis of TCD measurements. Combining these numbers yields a typical Reynolds number range of 175 Re 300. For convenience, we adopt Re = 250 in the sequel, which, in terms of the chosen reference length-scale and kinematic viscosity, corresponds to a velocity scale of = 0.388 m/s, well within the quoted range found in literature. In the sequel we will also consider the dependence of the fluid dynamics in relation to variation in the Reynolds number, particularly when we consider pulsatile flow.
In order to visualize the flow and forces that develop in the aneurysm, a number of options is available. We will start by choosing a more qualitative set of methods, i.e., three-dimensional and two-dimensional views of the velocity and shear stress. Here, we show results obtained at a resolution of 128 ×64 ×128. A quantitative approach, showing the effect of grid refinement will be followed in the next subsection.
Numerically computed velocity streamlines for a steady flow at Re = 250 are presented in Figure 4. By properly selecting the initial condition for the streamline at the inflow on the left-hand side of the domain, we can achieve both streamlines that pass through the section with the aneurysm, without actually entering the aneurysm bulb, as well as streamlines that display the rather complex vertical pattern that appears within the aneurysm. The geometry is colored with the value of the local pressure. A smooth transition from a high pressure to a low pressure can be observed, driving the flow inside the geometry.
To get an alternative impression of the velocity field, in Figure 5 we present the velocity vector field in a few cross-sections near the aneurysm bulb. The cross-sections are taken in yz planes, largely perpendicular to the mean flow direction. The flow in more or less cylindrical tube sections of the local vessel system shows similarity to a parabolic profile, reminiscent of Hagen-Poiseuille flow. We notice that the dominant flow still follows what used to be the non-diseased tract (Figure 5(c),(d)). However, also within the aneurysm there is considerable dynamics, showing regions of forward and backward flow, as illustrated in the velocity vectors. In the middle cross-sections (Figure 5(a),(b)) the flow dynamics is nicely illustrated with flow entering (green arrows) as well as exiting (blue arrows) the aneurysm in a complex vortical sweep; the flow partially comes back from the aneurysm after having circulated in it. The flow structure inside the aneurysm also leads to higher residence times of red blood cells, and possibly a reduced quality of circulation. This is correlated with regions of lower levels of wall-shear stress, as the flow-intensity in the aneurysm is rather low. From this general visualization one may already appreciate the frequently expressed observation that regions of low wall-shear stress are a marker of increased likelihood of aneurysm growth, possibly since the endothelial cells lining the vessel walls are (slightly) less well supplied with oxygen and nutrients, leading to their gradual degeneration .
A closer impression of the flow inside the aneurysm can be obtained by considering contour plots of velocity components. In Figure 6 we show a particular contour of the streamwise velocity component at x = 8, which is through the actual aneurysm, i.e., a section through the middle in Figure 5(a). We show a cross-section in the yz-plane at a number of spatial resolutions. The grid refinement shows a clear qualitative convergence toward the grid-independent solution. The resolution 32 × 16 × 32 is seen to be insufficient to capture the full complexity of the flow. However, the main features are already captured properly at a resolution of 64×32×64, while high accuracy results can only be expected by further increase of the resolution. We notice both dark and light colors in this contour plot, corresponding to positive and negative streamwise velocities. These show a region of recirculating flow in the aneurysm, next to the main through-flow represented by the dark region in the lower left corner of each contour plot. We also investigated the dependence of the flow prediction on spatial resolution at other streamwise locations and found similar qualitative convergence. A more precise assessment of the level of convergence is considered momentarily.
Regions of relatively high and relatively low shear stress are considered as important markers for the risk of aneurysm growth. These can be obtained from the simulations as well, by post-processing the velocity field. The shear stress τ is defined in non-dimensional form as
where Sij = (∂iuj + ∂jui)/2 denotes the rate of strain tensor. A global impression of the wall shear stress distribution is given in Figure 7. Regions of high shear stress are concentrated near relatively sharp bends in the vessel and near the ‘neck’ of the aneurysm(red and yellow), where the bulge connects to the previously unaffected vessel. Inside the aneurysm the shear stress is rather low (uniform blue), consistent with the rather low velocities that are observed inside the aneurysm bulge. Such a region of low (wall) shear stress is reported to be connected to aneurysm growth, associated with a slow degeneration of endothelial cells at the vessel walls. Further research in this direction is highly needed to clarify the precise mechanisms and to quantify possible growth-paths of the aneurysm.
A more precise impression of the shear stress distribution can be obtained in terms of contour plots. In Figure 8 we show the steady state shear stress distribution in a yz-plane through the middle of the aneurysm at x = 8. We observe a qualitative convergence of the shear stress, although the degree of convergence seems to be slightly less compared to the velocity field as shown in Figure 6. For the shear stress we need to approximate the derivative of the velocity, which is more demanding on the spatial resolution, especially close to the aneurysm wall. The aneurysm region shows one focal point of somewhat higher shear stress, while elsewhere the level of the shear stress is seen to be rather low. In addition, quite high shear stress levels are observed in the lower left corner of each contour plot, corresponding to the main flow through what remains of the original vessel structure prior to the development of the aneurysm.
In order to quantify more precisely the level of convergence, in the next section we consider the reliability of the flow simulations by investigating the quality with which the actual local solution is obtained.
3.2. Reliability of IB predictions: a grid refinement study
In this subsection we concentrate on a more quantitative analysis of the reliability of the numerical flow prediction. We consider the convergence of the driving pressure drop as well as profiles of velocity and shear stress at different grid resolutions.
In order to maintain a constant mass flow rate through the aneurysm, a pressure drop ΔP needs to be supplied. In Figure 9 we show the evolution of the forcing pressure drop at different spatial resolutions. The initial solution from which each simulation starts uses u = 1 and v = w = p = 0. Hence, we deliberately set the streamwise velocity equal to unity everywhere, i.e., also inside the solid part of the domain. This presents a strict test for the robustness of the method, in which the solution in the solid has to adapt and reduce completely to zero within the solid. Moreover, a realistic flow in the fluid part needs to build up. There is considerable difference between the solution at different spatial resolutions in the initial stage due to the strong acceleration of the flow to rectify the non-physical aspects of the initial condition. As the flow settles into the steady state we notice a clear convergence of the pressure drop levels. Since the non-dimensional size of the domain is 16 and the velocity is maximally on the order of 0.7, a typical flow-through time, i.e., the time needed to pass from one side of the domain to the other, can be expected to be in the range of 20 or more. To reach a fully steady state, a simulation covering several flow-through times is needed. In Figure 9 we notice already a fair agreement for ΔP at different spatial resolutions after about one flow-through time.
To further assess the reliability of the solution we turn to profiles of velocity and shear stress in a characteristic region in Figure 10. This provides a quantitative measure for the convergence of the numerical solution. We observe that the coarsest resolution of 32×16×32 is not capable to capture more than the main flow pattern of recirculating flow. Increasing the resolution shows a much closer correspondence between the different velocity predictions. This is also seen in the profiles for the shear stress. The general agreement is quite close, as long as we do not include the coarsest resolution. Convergence of the sharp stress peak near the lower wall in this particular profile is seen to be most challenging to our IB method. We also investigated convergence by considering profiles in a range of other locations and observed similarly close agreement of the numerical solutions at different spatial resolutions. This establishes that a first quantitatively acceptable solution can be obtained using a grid of 64 × 32 × 64, while higher accuracy requires further refinement.
In the next section we turn our attention to realistic pulsatile flow in the given cerebral aneurysm geometry.
4. Pulsatile flow simulations in realistic cerebral aneurysm
In this section we will focus on pulsatile flow. We first consider the pulsatile velocity pattern and translate this into the volumetric flow rate used in the simulations. Then we will present the dynamics of the shear stress at different, physiologically relevant Reynolds numbers. Next to Re = 250 as used up to now, we include Re = 200 and Re = 300 to probe the variability of predictions when alternative values for the blood viscosity, the vessel sizes and/or velocity scales are taken from literature. We also compute the flow at Re = 100 and Re = 400, which are expected to be indicative of diseased conditions. A striking transition to very complex time-dependence at higher Re is observed, which may be of interest for medical monitoring.
The velocity of the blood flow in cerebral arteries can be measured by different techniques, e.g., phase-contrast (MR) angiography  or TCD sonography as presented in . In the current study the velocity was recorded in the middle cerebral artery by  using TCD. In Figure 11(a) we plotted a segment of 10 seconds of the pulsating velocity wave. The mean velocity value, obtained by integrating this signal, is found to be 38.66 cm/s, which is very close to the reference scale selected in Section 3.
The computed pulsatile flow is maintained by using the actually recorded velocity signal as forcing. We choose a typical pulse from Figure 11(a) and repeat it periodically. We need to convert the recorded velocity values into a time-dependent volumetric flow rate, which we specify next. The selected pulse has a maximal velocity which corresponds to a peak flow rate of using the selected radius R∗ = 1.94 mm and assuming a perfectly circular cross section. If we take the reference velocity corresponding to a Reynolds number Re = 250, we find similarly as reference flow rate For convenience, we split the forcing signal in the non-dimensional formulation into a normalized flow rate pattern Q0 and an amplitude such that the forcing used in the simulations becomes The normalized pulsatile wave Q0 is illustrated in Figure 11(b). The physical duration of one pulse is t∗ = 0.82 s. The reference time-scale can be computed as Thus at Re = 250 one pulse requires 164 non-dimensional time units.
The procedure to define the pulsatile flow rate can be extended to also address other Reynolds numbers. We take as reference Reynolds number Re and fix the reference length-scale to R∗ (since we consider the same geometry) and the kinematic viscosity to (since we still consider the flow of blood). If we wish to simulate at another Reynolds number Re this implies that the reference velocity scale is changed according to Correspondingly, the time-scale changes into and hence, the ‘new’ number of dimensionless time-steps to take in order to complete one cycle of 0.82 s of the pulsatile flow decreases with decreasing Reynolds number. Another consequence of changing the Reynolds number at constant length-scale and kinematic viscosity is that as well as Hence, the dimensionless forcing does not alter with changing Reynolds number and remains at Q(t) ≈ 1.75Q0(t). The factor 1.75 denotes the ‘contrast’ in the pulsatile flow rate, i.e., the ratio between the maximal and the average velocity during a cycle - this quantity varies from one person to another and even per heartbeat.
We compute the pulsatile flow in a range of Reynolds numbers 100 ≤ Re ≤ 400 while keeping the viscosity of the blood and the size of the vessel constant. In Section 3 we discussed the uncertainty in physical parameters, when consulting literature. This leads to a range of Reynolds number 175 Re 300 of physiologically realistic values. Including also unhealthy changes in the blood vessels, e.g., narrowing of the vessel diameter, or the development of an aneurysm and corresponding increase in the size of the vessel, but also changes in the viscosity of blood, or in the velocity of the flow, we propose to also compute the flow at Re = 100 and Re = 400. This total range of Reynolds numbers gives a complete set of flow conditions relevant to flow in the Circle of Willis. For all simulations we use a spatial resolution of 64 × 32 × 64, which we showed to be the lower limit at which quantitatively reliable results can be obtained for the case considered here.
The translation ‘back’ from non-dimensional units to physical units requires scaling of the time and of the shear stress values. For the indicated range of Reynolds numbers Re a single pulse of 0.82 s requires #t dimensionless time steps with (Re, #t) = (100, 65.6), (200, 131.2), (250, 164), (300, 196.8) and (400, 262.4). Moreover, the change in Re corresponds to a change in which affects the scale for the shear stress which is The final result is a wall shear stress in Pa and time measure in s, which allows a more direct assessment than the fully dimensionless representation.
After these preparations, we show the dynamic response of the maximum shear stress in Figure 12 for 4 full pulse cycles. The initial condition is that of quiescent flow, i.e., u = v = w = 0 - we observe that this leads to a short transient, e.g., seen because the response in the first cycle differs slightly from that in later cycles. After this transient we observe for our reference case Re = 250 (solid bold line) that the maximum shear stress closely follows the imposed volumetric flow rate profile. The mean value is found to be around 1.4 Pa with peak values near 2.9 Pa. These values show the same general magnitude as reported in [14, 37].
Increasing or decreasing the Reynolds number has a marked effect on the dynamic response. At the lower Reynolds numbers the response of the shear stress maximum is seen to be smooth, following the imposed pulsatile profile. At the higher Reynolds numbers the natural Navier-Stokes nonlinearity seems to become dominant, which makes the shear stress response lively by the emergence of relatively high frequency components to the solution. In addition, the level of the shear stress rises considerably. Such rapid transition in flow regime within the physiologically relevant Re-range may contribute to an increased risk of aneurysm rupture. This transition was also seen in simulations of other realistic aneurysms and even for simplified model aneurysms consisting of curved vessel to which a spherical cavity was added .
5. Concluding remarks
We presented computational modeling of cerebral flow based on a volume penalizing IB method, aimed at understanding the flow and forces that emerge in aneurysms that may form on the Circle of Willis. We sketched how medical imagery can serve as point of departure for a sequence of numerical representations and modeling steps, ultimately leading to the full simulation of pulsatile flow in a realistic cerebral aneurysm, which was used as a case study in this chapter.
Taking data from literature we identified physiologically relevant flow conditions and their general uncertainty. Data concerning sizes of vessels, kinematic viscosity of blood and flow speeds in the region of the Circle of Willis during the cardiac cycle, can not be obtained with very high accuracy. This leaves considerable uncertainty as to the precise flow conditions. However, consensus seems to exist that the Reynolds number, which is the crucial parameter for incompressible flow, should be in the range 175 ≤ Re ≤ 300 for non-diseased situations. This is a rather wide range, but throughout this range the flow in the blood vessels is pulsatile and laminar, i.e., sharing quite comparable dynamics. The main challenge for computational modeling is not just to predict a certain flow under specified conditions, but to reckon also with the variability in flow conditions due to a variety of possible changes in key parameters. This approach was taken in this chapter.
Settling for Re = 250 as characteristic point of reference, we analyzed a particular, realistic cerebral aneurysm in detail. First, the main flow features and the reliability of predictions were considered for steady flow at fixed volumetric flow rate. This is not a realistic flow condition, as in reality interest is with pulsatile flow, but it does allow to investigate the sensitivity of the predicted solution on things such as spatial resolution. We visualized both qualitatively and quantitatively the steady flow in the aneurysm, as well as the shear stress field that emerges. It was shown that the main flow follows a path that is close to what used to be the original vessel before the formation of the aneurysm. Next to this ‘main’ flow, a complex circulation was shown to develop inside the aneurysm bulge. By considering contour plots and also profiles of velocity and shear stress at different spatial resolutions, the degree of reliability of the numerical simulation was discussed. The current IB method is first order accurate. Developments in which sub-grid forcing is included  can be used to increase the formal order of accuracy to two - this appears a relevant extension of the IB approach and will be considered in more detail in the near future, allowing to cut down on the computational cost and/or increase the accuracy of flow predictions.
A complete model was obtained by incorporating realistic pulsatile flow, obtained from direct TCD measurements of the velocity of blood flow in the brain. The recorded velocity in a vessel near the Circle of Willis was used to impose a proper time-dependent volumetric flow rate, representing a cardiac cycle. Repeating a characteristic pulse periodically, leads to a model with which the time-dependent flow and shear stresses can be determined. As regards the flow structures, for pulsatile flow at Re = 250 one basically notices the same dominant flow features as in case of steady flow forcing, with the exception that the ‘amplitude’ of the motion becomes time-dependent. As an example, the recirculating flow in the aneurysm was observed throughout the entire cycle, but with a time-dependent intensity and a slight ‘meandering’ of the precise vortical structures. We sketched the extension of the pulsatile flow model to other flow conditions, i.e., other Re and other time-scales. If the flow is in the physiologically relevant regime Re ≤ 300, the response of, e.g., the shear stress, closely follows that of the input flow-rate forcing. At higher Reynolds numbers, indicative of possible diseased states, the flow develops considerable complexity and shows a transition toward much higher frequencies. This goes hand in hand with increased levels of shear stress and may be monitored as a potential indication of increased risk to the patient. By recording the spectrum of these frequencies an easy monitoring concept may become available.
The application of computational support in the monitoring and treatment of cerebral aneurysms is a field of ongoing research. Accessibility of time-dependent flow fields in all relevant detail is a crucial point from which to depart toward developing predictive capability for the associated slow growth of the aneurysm bulge. This requires a new kind of ‘flow-structure interaction’ in which degradation of endothelial cells due to reduced quality of blood circulation typically triggers a further expansion of the aneurysm bulge, and generally leads to an increase in the risk of rupture. The latter type of ‘flow-structure interaction’ is subject of ongoing research. In order to achieve a closer connection with medical practice several computational modeling steps still need to be taken, such as the development of higher order accurate methods, multi inflow/outflow configurations, flow-structure interaction in which a full coupling to slow degenerative processes of endothelial cells is made, as well as modeling of mechanical properties of brain tissue to also address aneurysm compliance during the pulsatile cycle. This brief listing shows the various developments that are still needed in order to make the pathway from medical imagery to quantitative decision support both reliable as well as fully automated.AppendixAcknowledgement
The authors gratefully acknowledge many fruitful discussions with Prof. Dr. Hans Kuerten (Eindhoven University of Technology and University of Twente). We are also grateful to Willem Jan van Rooij, MD, PhD and Menno Sluzewski, MD, PhD (St. Elizabeth Hospital, Tilburg, the Netherlands) for providing angiographic data and to Dr. Ir. Dirk-Jan Kroon (Focal, Oldenzaal, the Netherlands) for segmentation of the data. Computations at the Huygens computer at SARA were supported by NCF through the project SH061.