**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 [50]. 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 [40]. In one approach, the blood is approximated as a Newtonian fluid [5]. 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*m*. A length-scale that characterizes the cross section of a typical cerebral vessel inside the Circle of Willis is on the order of 10−3*m* [19]. This difference in length-scale of three orders of magnitude motivates to approximate blood as an incompressible Newtonian fluid [40].

4 Will-be-set-by-IN-TECH

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

For a more complete modeling of the dynamics in the vessel system, flow-structure interaction often plays a role [40]. 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 [35]. For realistic pulsatile flows some movement of the aneurysm walls was observed during a cardiac cycle [36]. 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

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

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 [50]. 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

There are various approaches to model flow of blood in the human brain. A comprehensive overview is given in [40]. In one approach, the blood is approximated as a Newtonian fluid [5]. 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

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*m*. A length-scale that characterizes

smoothing of the geometry or sophisticated reconstruction methods.

cerebral aneurysms for which the wall movement can be neglected [24].

presented in Section 4. Concluding remarks are in Section 5.

precise geometry of realistic blood vessels.

fluid at somewhat different stress and velocity levels.

**2. Numerical model of blood flow inside human brain**

**2.1. The Navier-Stokes equations for an incompressible cerebral flow**

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 Ω*<sup>f</sup>* 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:

$$\frac{\partial \mathbf{u}^\*}{\partial t^\*} + \mathbf{u}^\* \cdot \nabla^\* \mathbf{u}^\* = -\nabla^\*(\frac{p^\*}{\rho^\*}) + \nu^\* \nabla^{\*2} \mathbf{u}^\* + \frac{\mathbf{f}^\*}{\rho^\*} \tag{1}$$

$$\nabla^\* \cdot \mathbf{u}^\* = 0 \tag{2}$$

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 *u*∗ *<sup>r</sup>* and reference length *L*<sup>∗</sup> *<sup>r</sup>* we can express a reference time scale as *t*<sup>∗</sup> *<sup>r</sup>* = *L*<sup>∗</sup> *<sup>r</sup>* /*u*<sup>∗</sup> *r* . Using as reference density *ρ*∗ *<sup>r</sup>* = *ρ*<sup>∗</sup> we will use a reference pressure as *p*<sup>∗</sup> *<sup>r</sup>* = (*u*<sup>∗</sup> *<sup>r</sup>* )2*ρ*<sup>∗</sup> *<sup>r</sup>* . For the forcing term we select a direct volume penalization in which

$$\frac{\mathbf{f}^\*}{\rho^\*} = -\frac{1}{\varepsilon^\*} H \mathbf{u}^\* \tag{3}$$

where *ε*<sup>∗</sup> is a forcing time scale and *H* is the masking function: *H*(**x**) = 1 if **x** ∈ Ω*<sup>s</sup>* and *H*(**x**) = 0 if **x** ∈ Ω*<sup>f</sup>* . We set the reference forcing time scale *ε*<sup>∗</sup> = *εt* ∗ *<sup>r</sup>* � *t* ∗ *<sup>r</sup>* , 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:

$$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = -\nabla p + \frac{1}{Re} \nabla^2 \mathbf{u} - \frac{1}{\varepsilon} H \mathbf{u} \tag{4}$$

$$\nabla \cdot \mathbf{u} = 0 \tag{5}$$

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 *ε* [26]. 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 [53], which implies a smooth velocity and pressure field. With increasing Reynolds number the flow can develop more detailed vortical 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 [12]. 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 [12]. 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.

**Figure 1.** Sketch of a primary grid cell in 2D with staggered allocation of the variables. The pressure *p* is in the middle of the grid cell, while the velocities (*u* and *v*) are defined at the centers of the faces.

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

$$\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} - \nabla p + \frac{1}{Re} \nabla \cdot \nabla \mathbf{u} \tag{6}$$

we are interested in the kinetic energy, given by

$$E = \frac{1}{2} \int\_{\Omega} dV \left| \mathbf{u} \right|^2 = \frac{1}{2} \int\_{\Omega} dV \mathbf{u} \cdot \mathbf{u} \tag{7}$$

where **u** · **u** is the vector inner product. Note, that in (7) we effectively integrate only over Ω*<sup>f</sup>* as **u** = **0** in Ω*s*. The evolution of the kinetic energy follows from

$$\frac{dE}{dt} = \int\_{\Omega} dV \mathbf{u} \cdot \frac{\partial \mathbf{u}}{\partial t} \tag{8}$$

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

6 Will-be-set-by-IN-TECH

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

We employ a staggered allocation of the flow variables (**u**, *p*)=(*u*, *v*, *w*, *p*) as basis for our flow solver [12]. 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

vi,j

xi-1 x<sup>i</sup>

**Figure 1.** Sketch of a primary grid cell in 2D with staggered allocation of the variables. The pressure *p* is in the middle of the grid cell, while the velocities (*u* and *v*) are defined at the centers of the faces.

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.

*dE dt* <sup>=</sup> Ω *dV***u** ·

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>(**<sup>u</sup>** · ∇)**<sup>u</sup>** − ∇*<sup>p</sup>* <sup>+</sup>

<sup>2</sup> <sup>=</sup> <sup>1</sup> 2 Ω

where **u** · **u** is the vector inner product. Note, that in (7) we effectively integrate only over Ω*<sup>f</sup>*

1

*∂***u**

*Re* ∇·∇**<sup>u</sup>** (6)

*dV***u** · **u** (7)

*<sup>∂</sup><sup>t</sup>* (8)

Starting from the original momentum equation without the forcing term

*∂***u**

*<sup>E</sup>* <sup>=</sup> <sup>1</sup> 2 Ω *dV*|**u**|

as **u** = **0** in Ω*s*. The evolution of the kinetic energy follows from

we are interested in the kinetic energy, given by

pi,j

ui,j

**2.2. Numerical method for flow simulation using an immersed boundary**

**approach for representing complex geometries**

specify the volume penalization IB method afterwards.

as the velocity- and the pressure-points, respectively.

yj-1

yj

$$\frac{dE}{dt} = -\frac{1}{Re} \int\_{\Omega} dV (\nabla \mathbf{u} : \nabla \mathbf{u}) \le 0 \tag{9}$$

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 [50], which departs from the expressions

$$E = \frac{1}{2}(\mathbf{u}, \mathbf{u}) \quad \text{and} \quad \frac{dE}{dt} = \frac{1}{2}\frac{d}{dt}(\mathbf{u}, \mathbf{u}) \tag{10}$$

where energy is written in terms of the function inner product (**u**, **u**) as defined in (7). This yields the symmetric expression

$$\begin{split} \frac{dE}{dt} &= -\frac{1}{2} \Big( ( (\mathbf{u} \cdot \nabla) \mathbf{u}, \mathbf{u} ) + (\mathbf{u}, (\mathbf{u} \cdot \nabla) \mathbf{u}) \Big) - \frac{1}{2} \Big( (\nabla p, \mathbf{u}) + (\mathbf{u}, \nabla p) \Big) \\ &+ \frac{1}{2 \text{Re}} \Big( (\nabla \cdot \nabla \mathbf{u}, \mathbf{u}) + (\mathbf{u}, \nabla \cdot \nabla \mathbf{u}) \Big) \end{split} \tag{11}$$

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

$$\mathbf{A}\frac{d\mathbf{u}\_h}{dt} = -\mathbf{C}\mathbf{u}\_h - \mathbf{D}\mathbf{u}\_h + \mathbf{M}^T\mathbf{p}\_h\tag{12}$$

$$\mathbf{M}\mathbf{u}\_h = \mathbf{0} \tag{13}$$

where **u***<sup>h</sup>* is the vector containing the discrete velocity solutions (*uh*, *vh*, *wh*), **p***<sup>h</sup>* 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 <sup>−</sup>**M***T*, while the coefficient matrix **<sup>M</sup>** itself represents the discretization of the divergence operator, integrated over the control volumes [50].

The discrete approximation for the kinetic energy can be given using the midpoint rule, as

$$E\_h = \mathbf{u}\_h^T \mathbf{A} \mathbf{u}\_h \tag{14}$$

Similar to the continuous case (11) we compute the evolution of the energy in the discrete model as

$$\frac{d E\_h}{dt} = -\mathbf{u}\_h^T (\mathbf{C} + \mathbf{C}^T) \mathbf{u}\_h - \mathbf{u}\_h^T (\mathbf{D} + \mathbf{D}^T) \mathbf{u}\_h + \mathbf{u}\_h^T (\mathbf{M}^T \mathbf{p}\_h) + (\mathbf{M}^T \mathbf{p}\_h)^T \mathbf{u}\_h \tag{15}$$

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***<sup>T</sup>* = **0**. The two terms related to the numerical pressure gradient can be rewritten as

$$\mathbf{u}\_h^T(\mathbf{M}^T \mathbf{p}\_h) + (\mathbf{M}^T \mathbf{p}\_h)^T \mathbf{u}\_h = (\mathbf{M} \mathbf{u}\_h)^T \mathbf{p}\_h + \mathbf{p}\_h \mathbf{M} \mathbf{u}\_h = \mathbf{0} \tag{16}$$

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:

$$\frac{dE\_h}{dt} = -\mathbf{u}\_h^T (\mathbf{D} + \mathbf{D}^T) \mathbf{u}\_h \le 0 \tag{17}$$

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 [50]. 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 [43] 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 [50].

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 Ω*<sup>f</sup>* is not too complex and does not contain too many separate objects around which the fluid should flow [22]. 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 [27]. 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 [28].

We use an IB method based on volume penalization [32] 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−<sup>10</sup> in the sequel. As a result, the momentum equation (4) reduces to *∂t***u** ≈ −**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.

8 Will-be-set-by-IN-TECH

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

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 [50]. 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 [43] 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

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 Ω*<sup>f</sup>* is not too complex and does not contain too many separate objects around which the fluid should flow [22]. 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 [27]. 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 [28]. We use an IB method based on volume penalization [32] to capture flow in complex aneurysm geometries. A key role in our IB method is played by the so-called 'masking function' *H*,

*<sup>h</sup>* (**<sup>D</sup>** <sup>+</sup> **<sup>D</sup>***T*)**u***<sup>h</sup>* <sup>≤</sup> <sup>0</sup> (17)

*dEh*

*dt* <sup>=</sup> <sup>−</sup>**u***<sup>T</sup>*

well-known Adams-Bashforth scheme. More details can be found in [50].

provide a strict decrease of the energy:

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 [31], 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.

**Figure 2.** IB method illustrated on a Cartesian grid in 2D for an arbitrary geometry. Dark cells correspond to fluid cells with masking function *H* = 0, while white cells represent solid part of the computational domain, where *H* = 1.

#### **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 [23]. 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.

**Figure 3.** Three dimensional geometries of realistic aneurysms reconstructed from 3DRA. The segmented data are presented in (a), while in (b) we plot the 'cut'-vessel. The original aneurysm region (red) and the 'connecting' artificial vessel (black) are displayed.

#### **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.

10 Will-be-set-by-IN-TECH

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 [23]. 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

(a) (b)

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

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

**Figure 3.** Three dimensional geometries of realistic aneurysms reconstructed from 3DRA. The segmented data are presented in (a), while in (b) we plot the 'cut'-vessel. The original aneurysm region

and the reliability of the predictions quantified.

(red) and the 'connecting' artificial vessel (black) are displayed.

shear stress in Subsection 3.2.

**3. Steady flow simulations in realistic cerebral aneurysm**

**3.1. Motivation and exact definition of the reference case**

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 [11] who adopt 2.5 *mm*, and also consistent with [19], 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 *<sup>ρ</sup>*<sup>∗</sup> is in the range 1025 <sup>≤</sup> *<sup>ρ</sup>*<sup>∗</sup> <sup>≤</sup> <sup>1125</sup> *kg*/*m*3, while the dynamic viscosity of blood is reported to be 3 · <sup>10</sup>−<sup>6</sup> <sup>≤</sup> *<sup>μ</sup>*<sup>∗</sup> <sup>≤</sup> <sup>4</sup> · <sup>10</sup>−<sup>6</sup> *Pa s*. Specifically, choosing typical values for the mass density of blood *ρ*<sup>∗</sup> = 1060 *kg*/*m*<sup>3</sup> and the dynamic viscosity of blood *<sup>μ</sup>*<sup>∗</sup> <sup>=</sup> 3.2 · <sup>10</sup>−<sup>3</sup> *Pa s* implies a kinematic viscosity *<sup>ν</sup>*<sup>∗</sup> <sup>=</sup> *<sup>μ</sup>*∗/*ρ*<sup>∗</sup> <sup>=</sup> 3.01 · <sup>10</sup>−<sup>6</sup> *<sup>m</sup>*2/*s*. Finally, the reference flow-rate as proposed by [15] and [34] is 245 <sup>±</sup> <sup>65</sup> *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 *<sup>u</sup>*<sup>∗</sup> <sup>=</sup> *<sup>Q</sup>*∗/(*π*(*R*∗)2) = 0.345 <sup>±</sup> 0.09 *<sup>m</sup>*/*s*. This is consistent with the range 0.34 ± 0.087 *m*/*s* as obtained by [41] 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 *u*∗ *<sup>r</sup>* = 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

**Figure 4.** Three characteristic velocity streamlines inside the aneurysm geometry. The geometry is colored with the pressure field which shows a smooth transition from a high pressure (red) to a low pressure (blue) area. Grid resolution is 128 × 64 × 128.

display the rather complex vortical 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.

**Figure 5.** Velocity vector field in a few cross-sections along the aneurysm geometry. The cross-sections are taken at several *x* locations in *yz* planes, largely perpendicular to the flow direction. Green and yellow arrows show the middle range of flow velocities. Areas of slightly higher velocity are shown in red, while blue arrows indicate negative velocities related to the recirculating flow that develops. Grid resolution is 128 × 64 × 128.

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 [8].

12 Will-be-set-by-IN-TECH

**Figure 4.** Three characteristic velocity streamlines inside the aneurysm geometry. The geometry is colored with the pressure field which shows a smooth transition from a high pressure (red) to a low

display the rather complex vortical 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

(a) (b)

(c) (d)

**Figure 5.** Velocity vector field in a few cross-sections along the aneurysm geometry. The cross-sections are taken at several *x* locations in *yz* planes, largely perpendicular to the flow direction. Green and yellow arrows show the middle range of flow velocities. Areas of slightly higher velocity are shown in red, while blue arrows indicate negative velocities related to the recirculating flow that develops. Grid

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

pressure (blue) area. Grid resolution is 128 × 64 × 128.

resolution is 128 × 64 × 128.

low pressure can be observed, driving the flow inside the geometry.

**Figure 6.** Contour plot of the streamwise *u*-velocity in a *yz* cross-section in the middle of the geometry at *x* = 8. Dark regions correspond to high positive velocities, while the lightest contours are related to regions of negative velocity; we adopt the same color-coding for all figures. The flow structures show a vortex inside the aneurysm. We present the velocity contours in the same plane for different grid resolutions: 32 × 16 × 32 (a), 64 × 32 × 64 (b) and 128 × 64 × 128 (c).

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.

**Figure 7.** Distribution of the non-dimensional shear stress, computed at *Re* = 250. The grid resolution is 128 × 64 × 128. Red and yellow areas correspond to the locations of high shear stress, while blue areas are related to the low stresses.

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

$$
\pi = \frac{1}{\text{Re}} \sqrt{2S\_{ij}S\_{ij}} \tag{18}
$$

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.

**Figure 8.** Contour plot of shear stress *τ* in a *yz* cross-section in the middle of the geometry *x* = 8. Dark regions correspond to high levels of shear stress, while the lightest contours are related to low shear regions. We present the distribution of stress contours in the same plane for different grid resolutions: 32 × 16 × 32 (a), 64 × 32 × 64 (b) and 128 × 64 × 128 (c).

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**

14 Will-be-set-by-IN-TECH

**Figure 7.** Distribution of the non-dimensional shear stress, computed at *Re* = 250. The grid resolution is 128 × 64 × 128. Red and yellow areas correspond to the locations of high shear stress, while blue areas

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

y

(b)

**Figure 8.** Contour plot of shear stress *τ* in a *yz* cross-section in the middle of the geometry *x* = 8. Dark regions correspond to high levels of shear stress, while the lightest contours are related to low shear regions. We present the distribution of stress contours in the same plane for different grid resolutions:

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

<sup>5</sup> 4.5 <sup>4</sup> 3.5 <sup>3</sup> 2.5 <sup>2</sup> 1.5 <sup>1</sup> <sup>5</sup>

2*SijSij* (18)

z

5.5 6 6.5 7 7.5 8 8.5 9 9.5

y

(c)

<sup>5</sup> 4.5 <sup>4</sup> 3.5 <sup>3</sup> 2.5 <sup>2</sup> 1.5 <sup>1</sup> <sup>5</sup>

*<sup>τ</sup>* <sup>=</sup> <sup>1</sup> *Re* 

and to quantify possible growth-paths of the aneurysm.

z

5.5 6 6.5 7 7.5 8 8.5 9 9.5

y

32 × 16 × 32 (a), 64 × 32 × 64 (b) and 128 × 64 × 128 (c).

(a)

<sup>5</sup> 4.5 <sup>4</sup> 3.5 <sup>3</sup> 2.5 <sup>2</sup> 1.5 <sup>1</sup> <sup>5</sup>

z

5.5 6 6.5 7 7.5 8 8.5 9 9.5

are related to the low stresses.

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.

**Figure 9.** Convergence of the pressure drop across the whole flow domain, computed at different grid resolutions: 32 × 16 × 32 (dash-dot), 64 × 32 × 64 (dash) and 128 × 64 × 128 (solid).

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.

**Figure 10.** Streamwise velocity (a) and shear stress (b) profiles as a function of *z*, in the middle of the domain at *x* = 8 and *y* = 4. The profiles were computed at *t* = 20 using different grid resolutions: 32 × 16 × 32 (dash-dot), 64 × 32 × 64 (dash) and 128 × 64 × 128 (solid).

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 [15] or TCD sonography as presented in [41]. In the current study the velocity was recorded in the middle cerebral artery by [9] 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.

16 Will-be-set-by-IN-TECH

5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10

**Figure 10.** Streamwise velocity (a) and shear stress (b) profiles as a function of *z*, in the middle of the domain at *x* = 8 and *y* = 4. The profiles were computed at *t* = 20 using different grid resolutions:

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

In the next section we turn our attention to realistic pulsatile flow in the given cerebral

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 [15] or TCD sonography as presented in [41]. In the current study the velocity was recorded in the middle cerebral artery by [9] using TCD. In Figure 11(a) we plotted a segment of 10 seconds of the pulsating velocity wave. The mean

**4. Pulsatile flow simulations in realistic cerebral aneurysm**

z

<sup>0</sup> 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 4.5

τ

(b)

−0.2 −0.1 <sup>0</sup> 0.1 0.2 0.3 <sup>5</sup>

(a)

higher accuracy requires further refinement.

aneurysm geometry.

U

32 × 16 × 32 (dash-dot), 64 × 32 × 64 (dash) and 128 × 64 × 128 (solid).

5.5 6 6.5 7 7.5 8 8.5

z

**Figure 11.** Pulsatile wave measured in the human brain using TCD sonography. The 10 seconds velocity signal measured in the MCA is presented in (a). The unit pulse mass-flow wave was chosen from the original signal and normalized for the computations (b).

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 *u*∗ *max* ≈ 67.94 *cm*/*s*, which corresponds to a peak flow rate of *Q*∗ *max* <sup>≈</sup> 8.033 · <sup>10</sup>−<sup>6</sup> *<sup>m</sup>*3/*s*, using the selected radius *<sup>R</sup>*<sup>∗</sup> <sup>=</sup> 1.94 *mm* and assuming a perfectly circular cross section. If we take the reference velocity *u*∗ *<sup>r</sup>* = 0.388 *m*/*s* corresponding to a Reynolds number *Re* = 250, we find similarly as reference flow rate *Q*∗ *<sup>r</sup>* <sup>≈</sup> 4.59 · <sup>10</sup>−<sup>6</sup> *<sup>m</sup>*3/*s*. For convenience, we split the forcing signal in the non-dimensional formulation into a normalized flow rate pattern *Q*<sup>0</sup> and an amplitude *Q*∗ *max*/*Q*<sup>∗</sup> *<sup>r</sup>* such that the forcing used in the simulations becomes *Q*(*t*)=(*Q*∗ *max*/*Q*<sup>∗</sup> *<sup>r</sup>* )*Q*0(*t*) ≈ 1.75*Q*0(*t*). The normalized pulsatile wave *Q*<sup>0</sup> 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 *t* ∗ *<sup>r</sup>* = *R*∗/*u*<sup>∗</sup> *<sup>r</sup>* = 0.005 *s*. 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 *ν*∗ *<sup>r</sup>* (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 (*u*∗)� *<sup>r</sup>* = (*Re*� /*Re*)*u*∗ *r* . Correspondingly, the time-scale changes into (*t*∗)� *<sup>r</sup>* = (*Re*/*Re*� )*t*∗ *<sup>r</sup>* 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 (*Q*∗)� *<sup>r</sup>* = (*Re*� /*Re*)*Q*∗ *<sup>r</sup>* , as well as (*Q*∗)� *max* = (*Re*� /*Re*)*Q*∗ *max*. Hence, the dimensionless forcing does not alter with changing Reynolds number and remains at *Q*(*t*) ≈ 1.75*Q*0(*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 #*<sup>t</sup>* 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 *u*∗ *<sup>r</sup>* , which affects the scale for the shear stress which is *ρ*<sup>∗</sup> *<sup>r</sup>* (*u*<sup>∗</sup> *<sup>r</sup>* )2. 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.

**Figure 12.** Maximal shear stress for realistic pulsatile flow in the aneurysm geometry. The reference *Re* = 250 is shown (bold, solid). Lower values are shown as *Re* = 200 (dot), *Re* = 100 (dash-dot), while higher value are displayed as *Re* = 300 (dash) and *Re* = 400 (thin, solid). The average stress levels decrease with decreasing Reynolds number.

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 [30].
