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−6
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 Ω
Here u∗ is the velocity of the fluid,
where ε∗ is a forcing time scale and
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
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,
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 Ω
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 =
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
In a discrete setting the Navier-Stokes equations in matrix-vector notation are written as
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 + C
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 Ω
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’
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
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
The raw data of the 3DRA scan of the aneurysm consists of a grid of 2563 voxels. The voxel width is 0.1213
Next to the length-scale, also the viscosity and the velocity scales need to be set. From literature we infer that the mass density
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
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
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
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
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
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 Δ
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
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
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
The procedure to define the pulsatile flow rate can be extended to also address other Reynolds numbers. We take as reference Reynolds number
We compute the pulsatile flow in a range of Reynolds numbers 100 ≤
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
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.,
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
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 ≤
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
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.
Angot P. Bruneau C. H. Frabrie P. 1999 A penalization method to take into account obstacles in viscous flows. Numerical Mathematics, 81 497 520
Bernsdorf J. Wang D. 2009 Non-Newtonian blood flow simulation in cerebral aneurysms, 58 1024 1029
Bescós J. O. Slob M. J. Slump C. H. Sluzewski M. van Rooij W. J. 2005 Volume Measurement of Intracranial Aneurysms from 3D Rotational Angiography: Improvement of Accuracy by Gradient Edge Detection.American Journal of Neuroradiology, 26 2569 2572
Boussel L. Rayz V. Mc Culloch C. Martin A. Acevedo-Bolton G. Lawton M. Higashida R. Smith W. S. Young W. L. Saloner D. 2008 Aneurysm Growth Occurs at Region of Low Wall Shear Stress: Patient-Specific Correlation of Hemodynamics and Growth in a Longitudinal Study 39 2997 3002
Cebral J. R. Castro M. A. Burgess J. E. Pergolizzi R. S. Sheridan M. J. Putman C. M. 2005 Characterization of Cerebral Aneurysms for Assessing Risk of Rupture By Using Patient-Specific Computational Hemodynamics Models.American Journal of Neuroradiology, 26 2550 2559
Cebral J. R. Castro M. A. Appanaboyina S. Putman C. M. Millan D. Frangi A. F. 2005bEfficient Pipeline for Image-Based Patient-Specific Analysis of Cerebral Aneurysms Hemodynamics: Technique and Sensitivity. Imaging, 24 4 457 467
Cheny Y. Botella O. 2010 The LS-STAG method: A new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties 229 4 1043 1076
Doenitz C. Schebesch K. M. Zoephel R. Brawanski A. 2010 A mechanism for the rapid development of intracranial aneurysms: a case study 5 1213 1221
Duijvenboden S. Schaafsma A. Geurts B. J. 2011Usefulness of Schrodinger´s operator for assessing new parameters in Transcranial Doppler sonography, University of Twente, Internal Report
Gao T. Tseng Y. H. Lu X. Y. 2007 An improved hybrid Cartesian/immersed boundary method for fluid-solid flows 55 12 1189 1211
Gambaruto A. M. Janela J. Moura A. Sequeira A. 2011 Sensitivity of Hemodynamics in a Patient Specific Cerebral Aneurysm to Vascular Geometry and Blood Rheology 8 2 409 423 doi:10.3934/mbe.2011.8.409
Geurts B. J. 2003 Elements of Direct and Large-Eddy SimulationEdwards Publishing, 1-93021-707-2
Gijsen F. J. H. van de Vosse F. N. Janssen J. D. 1999 The influence of the non-Newtonian properties of blood on the flow in large arteries: steady flow in a carotid bifurcation model., 32 601 608
Goubergrits L. Schaller J. Kertzscher U. van der Bruck N. Poethkow K. Petz Ch. Hege H. Ch Spuler A. 2012 Statistical wall shear stress maps of ruptured and unruptured middle cerebral artery aneurysms 9 69 677 688 doi:rsif.2011.0490
Hendrikse J. van Raamt A. F. van der Graaf Y. Mali W. P. T. M. van der Grond J. 2005 Distribution of Cerebral Blood Flow in the Circle of Willis. 235 184 189 doi:radiol.2351031799
Hirsch C. 1988 Numerical computation of internal and external flowsWiley and Sons, 0-47191-762-1 0471 917621
Iaccarino G. Verzicco R. 2003 Immersed boundary technique for turbulent flow simulationsASME, 56 3 331 347 doi:
Janela J. Moura A. Sequeira A. 2010 A 3D non-Newtonian fluid-structure Interaction model for blood flow in arteries, 234 9 2783 2791 doi:j.cam.2010.01.032
Kamath S. 1981 Observations on the length and diameter of vessels forming the Circle of willis., 133 3 419 423
Keetels G. H. Órtona D. U. Kramer W. Clercx H. J. H. Schneider K. van Heijst G. J. F. 2007 Fourier spectral and wavelet solvers for the incompressible Navier-Stokes equations with volume-penalization: Convergence of a dipole-wall collision 227 2 919 945 doi:j.jcp.2007.07.036
Khadra K. Angot P. Parneix S. Caltagirone J. P. 2000 Fictitious domain approach for numerical modelling of Navier-Stokes equations 34 8 651 684 doi:
Kovalev K. 2005Unstructured hexahedral non-conformal mesh generation. PhD Thesis, Vrije University of Brussels
Kroon D. J. Slump C. H. Sluzewski M. van Rooij W. J. 2006Image based hemodynamics modelling of cerebral aneurysms and the determination of the risk of rupture. Images, 6143Manduca, A; Amini, A.A. (Ed.), Proc. of SPIE, 61432B. doi:
Ku D. N. 1997 Blood flow in arteries, 29 399 434
Koffijberg H. Buskens E. Algra A. Wermer M. J. H. Rinkel G. J. E. 2008 Growth rates of intracranial aneurysms: exploring constancyJournal of Neurosurgery, 109 2 176 185
Liu Q. Vasilyev O. V. 2007 A Brinkman penalization method for compressible flows in complex geometries 227 2 946 966 doi:j.jcp.2007.07.037
Löhner R. Baum J. D. Mestreau E. L. Rice D. 2007 Comparison of body-fitted, embedded and immersed 3-d euler predictions for blast loads on columns. American Institute of Aeronautics and Astrounatics, 1 22
Lopez Penha. D. J. Geurts B. J. Stolz S. Nordlund M. 2011 Computing the apparent permeability of an array of staggered square rods using volume-penalization& Fluids, 51 1 157 173 doi:j.compfluid.2011.08.011
Metcalfe R. W. 2003 The Promise of Computational Fluid Dynamics As a Tool for Delineating Therapeutic Options in the Treatment of Aneurysms.Editorials, 24 553 554
Mikhal J. Geurts B. J. 2011 Pulsatile flow in model cerebral aneurysms 4 811 820 doi:j.procs.2011.04.086
Mikhal J. Geurts B. J. 2011 Bounding Solutions for Cerebral Aneurysms, 5 3 163 168
Mittal R. Iaccarino G. 2005 Immersed Boundary Methods 37 239 261 doi:annurev.fluid.37.061903.175743
Moret J. Kemkers R. Op de Beek. J. Koppe R. Klotz E. Grass M. 1998 3D rotational angiography: Clinical value in endovascular treatment 42 3 8 14
Oktar S. O. Yücel C. Karaosmanoglu D. Akkan K. Ozdemir H. Tokgoz N. Tali T. 2006 Blood-flow volume quantification in internal carotid and vertebral arteries: comparison of 3 different ultrasound techniques with phase-contrast MR imaging.American Journal of Neuroradiology, 27 2 363 369
Oubel E. De Craene M. Putman C. M. Cebral J. R. Frangi A. F. 2007 Analysis of intracranial aneurysm wall motion and its effects on hemodynamic patterns. SPIE Medical Imaging 2007: Physiology, Function, and Structure from Medical Images, 6511Manduca, A; Hu, X.P. (Ed.), San Diego, CA: SPIE Press, 10-11
Oubel E. Cebral J. R. De Craene M. Blanc R. Blasco J. Macho J. Putman C. M. Frangi A. F. 2010 Wall motion estimation in intracranial aneurysms 31 1119 1135 doi:
Oyre S. Ringgaard S. Kozerke S. Paaske W. P. Erlandsen M. Boesiger P. Pedersen E. M. 1998Accurate Noninvasive Quantification of Blood Flow, Cross-Sectional Lumen Vessel Area and Wall Shear Stress by Three-Dimensional Paraboloid Modeling of Magnetic Resonance Imaging Velocity Data. Cardiology, 32 1 128 134
Perktold K. Peter R. Resch M. 1989 Pulsatile non-Newtonian blood flow simulation through a bifurcation with an aneurysm. 26 6 1011 1030
Peskin C. S. 2002 The immersed boundary method 1 39 doi:S0962492902000077
Quarteroni A. Formaggia L. 2004 Mathematical Modelling and Numerical Simulation of the Cardiovascular System, North-Holland, Amsterdam
Ringelstein E. B. Kahlscheuer B. Niggemeyer E. Otis S. M. 1990 Transcranial Doppler Sonography: Anatomical landmarks and normal velocity values. 16 8 745 761
Sarthou A. Vincent S. Caltagirone J. P. Angot P. . 2008 Eulerian-Lagrangian grid coupling and penalty methods for the simulation of multiphase flows interacting with complex objects 56 1093 1099 doi:fld.1661
Saylor P. E. 1988 Leapfrog variants of iterative methods for linear algebraic equations, 24 169 193
Shojima M. Oshima M. Takagi K. Torii R. Hayakawa M. Katada K. Morita A. Kirino T. 2004 Magnitude and Role of Wall Shear Stress on Cerebral Aneurysm. Computational Fluid Dynamic Study of 20 Middle Cerebral Artery Aneurysms., 35 2500 2505 doi:STR.0000144648.89172.0f
Sprengers M. E. Schaafsma J. van Rooij W. J. Sluzewski M. Rinkel G. J. E. Velthuis B. K. van Rijn J. C. Majoie C. B. 2008 Stability of Intracranial Aneurysms Adequately Occluded 6 Months after Coiling: a 3T MR Angiography Multicenter Long-Term Follow-Up Study, 29 1768 1774 doi:ajnr.A1181
van Rooij W. J. 1998 Endovascular Treatment of Cerebral AneurysmsPhD Thesis, University of Utrecht. 9-03931-818-2
van Rooij W. J. Sprengers M. E. Sluzewski M. Beute G. N. 2007 Intracranial aneurysms that repeatedly reopen over time after coiling: imaging characteristics and treatment outcome. 49 343 349 doi:s00234-006-0200-2
van Rooij W. J. Sprengers M. E. de Gast A. N. Peluso A. P. P. Sluzewski M. 2008 3D Rotational Angiography: The New Gold Standard in the Detection of Additional Intracranial AneurysmsAJNR 29 976 979 doi:ajnr.A0964
Vasilyev O. V. Kevlahan N. K. R. 2002 Hybrid wavelet collocation-Brinkman penalization method for complex geometry flows 40 3-4 531 538 doi:fld.307
Verstappen R. W. C. P. Veldman A. E. P. 2003Symmetry-Preserving Discretization of Turbulent Flow. Journal ofComputational Physics, 187 343 368
Wanke I. Dörfler A. Forsting M. 2008 Intracranial Aneurysms. Intracranial Vascular Malformations and Aneurysms. Springer, 978-3-54032-919-0
Wiebers D. O. Whisnant J. P. Forbes G. et al. 1998Unruptured intracranial aneurysms- risk of rupture and risks of surgical intervention. International Study of Unruptured Intracranial Aneurysms Investigators. , 339 24 1725 1733
Young D. F. Munson B. R. Okiishi T. H. 1997 A brief introduction to fluid mechanicsJohn Wiley & Sons, Inc., 0-47113-771-5