**4. Modeling of blood flow around intracranial aneurysm**

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

**Figure 5.** Examples of our simulation results: (left) real coil embolization; (right) our simulated coil

*d***F** *hd***<sup>v</sup>** <sup>+</sup>

**<sup>W</sup>**(*α*,*β*) <sup>=</sup> *<sup>g</sup>*(*dij*)

*d***F** *d***q**

where **C**(*i*,*j*) is the 3 × 3 sub-matrix of global compliance matrix **C** (inverse of tangent matrix) at the rows of node *i* and the columns of node *j*. For the aneurysm wall, the formulation of

*e*

where *e* is an elasticity parameter that is homogeneous to young modulus and *g*(*dij*) is a Gaussian function of the distance, defined on the surface, between contact point *i* and *j*. The Gaussian function allows a fall-off of the coupling with increasing distance between the contact points. This model is close to the Boussinesq approximation which provides a distribution of the normal contact stress from the elasticity of the surface, around a point of

The result of the *contact response* involves finding the friction contact forces that respect Signorini and Coulomb laws. Several works ([18] and [8]) present Gauss-Seidel iterative approaches that solve this problem. The solver needs an evaluation of a global compliance matrix **W**, which is the sum of the compliances of the coil and the aneurysm wall. It also needs the value of the relative displacement of the contacting points during the free motion δfree. When the contact forces are found, during the last step, called *contact correction* we compute

In this section, we presented an efficient dynamic coil model, *i.e.*, FEM beams based model for the intrinsic mechanical behavior of the coil, combined with a process for computing collision response (with the aneurysm wall) and auto-collision response. However, even if the results are quite encouraging, as illustrated in Figure 5, the mechanical modeling is not fully complete. Indeed, the behavior of the coil is also greatly influenced by the blood flow.

−<sup>1</sup> (*i*,*j*)

**H***<sup>T</sup> <sup>α</sup>* **H***β*,

**<sup>H</sup>***<sup>β</sup>* = **<sup>H</sup>***<sup>T</sup>*

*<sup>α</sup>* **C**(*i*,*j*)**H***β*,

surface and the points of self-collision. When collisions are detected, the *contact response* is computed. This is a complex aspect that influences greatly the overall behavior of the coil model. To describe the mechanical behavior during contact, the mechanical coupling between the different contact points is modeled. This information is provided by evaluating the compliance matrix in the contact space, called **W**, for both the coil and the aneurysm. Let's consider a contact *α* on the node *i* of the coil (with one constraint along the contact normal **n** and two along the tangential friction directions **t**, **s**). **H***<sup>α</sup>* is the matrix of the frame [**nts**]. The mechanical coupling of this contact with a contact *β* (with frame **H***β*) on node *j* can be

embolization with 3D coils.

the coupling is simpler:

contact [30].

evaluated with the following 3 × 3 matrix:

the motion associated to the contact forces.

**<sup>W</sup>**(*α*,*β*) <sup>=</sup> **<sup>H</sup>***<sup>T</sup>*

*α* **M** *<sup>h</sup>*<sup>2</sup> <sup>+</sup> The blood is modeled as an incompressible Newtonian fluid of constant density, *ρ* = 1.069 × 103*kg*/*m*3, and constant viscosity, *<sup>μ</sup>* <sup>=</sup> 3.5 <sup>×</sup> <sup>10</sup>−3*kg*/(*<sup>m</sup>* ·*s*). In this work, we only consider the blood flow near intracranial aneurysms with relatively small Reynolds number ranging from 100 to 1,000 [31], which satisfies the laminar assumption. Thus, it is reasonable to describe the behavior of blood flow using the incompressible Navier-Stokes equations. Additionally, the influence of pulsating vessel on the fluid is ignored; hence vessel walls are assumed to be rigid.

We rely on the Discrete Exterior Calculus (DEC) method for numerically solving the equations, as it offers several benefits for our application. First of all, it is based on unstructured tetrahedral mesh, which is more suitable to describe irregular boundary of anatomical geometries than regular grids. From the tetrahedral mesh, referred as primal mesh, the dual mesh is constructed as follows: dual vertices correspond to the circumcenters of primal tetrahedra, dual edges link dual vertices located on neighbor tetrahedra, and dual faces are surfaces of Voronoi cells of primal vertices, which are dual cells as well. More generally, a dual (*n* − *p*)-cell is associated to a corresponding *p*-simplex (*p* = 0, 1, 2, 3, *n* = 3) as depicted in Figure 7.

Secondly, the orthogonality between primal and circumcentric dual is useful to define the physical variables in fluid (Figure 6), such as flux, which is the face integral of velocity orthogonal to the face, and the vorticity (the circulation per unit area at a point), whose direction is orthogonal to the plane of circulation (the line integral of velocity around a closed curve). The discrete version of these physical quantities is not only defined on points, but represented as scalars attached to the primitives of any dimension on primal or dual mesh, *i.e.*, vertex, edge, face or volume. The scalars are integral values of physical quantity over the primitives (except point-wise variable), and the orientation of each primitive represents the local direction of vector variable. For example, velocity is described as flux on each triangle face (2D primal primitive), so it is primal 2-form, and the orientation of the face indicates whether the fluid flows outward or inward. Vorticity is defined as the integral value over the dual face (2D dual primitive), thus it is dual 2-form. According to the Stoke's theorem, this value equals to the circulation along the boundary of dual face. The orientation of the primal edge gives the direction of vorticity.

Thirdly, based on the discretization of space and variables, discrete vector calculus operators, such as curl, divergence, Laplace operators, can be easily expressed by two basic operators, the discrete differentials *d* and the Hodge stars *�*. The former, *dp*, maps *p*-forms to (*p* + 1)-forms on the primal mesh, represented by signed incidence matrix. The transpose of this matrix operates similarly on the dual mesh. The latter, *�p*, maps from primal *p*-forms to dual (*n* − *p*)-forms, represented by a diagonal matrix whose element equals to the volume ratio between the corresponding dual and primal elements. And the inverse matrix maps in the opposite (Figure 7).

**Figure 7.** The duality of primal and dual mesh: the first line shows the primal simplex, whose dual elements are below. Physical variables **U** and **Ω**, defined as discrete forms, can be transferred by two fundamental operators *d* and *�*.

Finally, it applies the same idea of Stable Fluids [39], a semi-Lagrangian method, which is much more stable than traditional Eulerian methods, and allows larger time steps and improves the computational efficiency. But instead of dealing with velocity field, the DEC method uses vorticity-based formulations (Equation 2), and preserves the circulation at a discrete level. As a result, to some degree, it overcomes the issue of numerical diffusion, and results in higher accuracy.

$$\begin{aligned} \frac{\partial \omega}{\partial t} &= -\mathcal{A}\_v \omega + \frac{\mu}{\rho} \Delta \omega \\ \nabla v &= 0 \quad \omega = \nabla \times v, \end{aligned} \tag{2}$$

where ω is the vorticity, Lie derivative Lvω (in this case equal to v · ∇ω − ω · ∇v), is the advection term, and *<sup>μ</sup> <sup>ρ</sup>* Δω is the diffusion term.

Simply speaking, the advection term describes the idea that the local spin is pushed forward along the direction of the velocity, which is consistent with Kelvin's circulation theorem: the circulation around a closed curve moving with the fluid remains constant with time. In this

**Figure 8.** Backtracking. At the current time step *tn* = *t*, given the velocity field, each dual vertex (on the right) is backtracked to its position at the previous time step *tn*−<sup>1</sup> = *t* − *h* (on the left). The circulation around the loop of dual face boundary at time *tn* is forced to be equal to the circulation of the backtracked loop at time *tn*−1, *i.e.*, *Ω*(*t*) = *Ω*(*t* − *h*).

approach, the discrete vorticity is conserved by extending Kelvin's theorem to the discrete level: the circulation around the loop of each dual face's boundary keeps constant as the loop is advected by fluid flow. So we run a backtracking step (Figure 8) to find out where the current dual face comes from, and accumulate the circulation around the backtracked dual face, and then assign this value to the current one. This step makes the computation circulation-preserving at a discrete level, as well as stable, because the maximum of the new field is never larger than that of the previous field. For the diffusion term, linear solver is used, and an implicit scheme is chosen for the purpose of stability (Equation 3).

$$
\omega\_{n+1} - \omega\_n = \frac{\mu \hbar}{\rho} L \omega\_{n+1} \tag{3}
$$

where *h* is the span of time step, and *L* is the Laplace operator.

#### **5. Blood-coil interaction**

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

(a) Velocity *v* and flux **U** (b) Vorticity *ω* and circulation **Ω**

**Figure 7.** The duality of primal and dual mesh: the first line shows the primal simplex, whose dual elements are below. Physical variables **U** and **Ω**, defined as discrete forms, can be transferred by two

*∂*ω

*<sup>ρ</sup>* Δω is the diffusion term.

Finally, it applies the same idea of Stable Fluids [39], a semi-Lagrangian method, which is much more stable than traditional Eulerian methods, and allows larger time steps and improves the computational efficiency. But instead of dealing with velocity field, the DEC method uses vorticity-based formulations (Equation 2), and preserves the circulation at a discrete level. As a result, to some degree, it overcomes the issue of numerical diffusion,

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>Lv<sup>ω</sup> <sup>+</sup> *<sup>μ</sup>*

∇v = 0 ω = ∇ × v,

where ω is the vorticity, Lie derivative Lvω (in this case equal to v · ∇ω − ω · ∇v), is the

Simply speaking, the advection term describes the idea that the local spin is pushed forward along the direction of the velocity, which is consistent with Kelvin's circulation theorem: the circulation around a closed curve moving with the fluid remains constant with time. In this

*ρ* Δω

(2)

fundamental operators *d* and *�*.

and results in higher accuracy.

advection term, and *<sup>μ</sup>*

**Figure 6.** Discretization of physical variables in fluid: (a) velocity; (b) vorticity.

In this section, we illustrate how we simulate the bilateral interaction between coils and blood flow. First, we describe how we compute the impact of the coil onto the blood flow by adding extra terms to the Navier-Stokes equation. Second, we explain how the drag force applied by blood flow onto the coil is computed. Finally, we combine the fluid and structure systems by a loosely-coupled strategy for real-time simulation of coil embolization.

#### **5.1. Porous media model**

The typical diameter of coils chosen for intracranial aneurysms ranges from 0.2*mm* to 0.4*mm*, which is quite small compared to the aneurysm size (0.02 to 0.11 times of the intracranial aneurysm in diameter). And after inserted into aneurysm, they are randomly distributed, forming the shape of a twisted nest (Figure 9). Considering the relatively small dimension of coils and their random distribution in aneurysm, coils are modeled, from a statistical point of view, as porous media in aneurysm.

We divide the fluid domain D into 3 sub-domains, a coil-free and a coil-filled sub-domain, as well as a transitory sub-domain between them, which allows the porous parameters to vary

**Figure 9.** The shape of detachable coils after deployment.

smoothly between the first two sub-domains. However, blood motion in all sub-domains is described uniformly by the Navier-Stokes equations (vorticity-based) of Brinkmann type:

$$\begin{aligned} \frac{\partial \boldsymbol{\omega}}{\partial t} + \mathcal{E}\_{\boldsymbol{\nu}} \boldsymbol{\omega} &= \frac{\mu}{\rho} \Delta \boldsymbol{\omega} - \frac{\mu \boldsymbol{\rho}}{\rho \boldsymbol{k}} \boldsymbol{\omega} - \frac{\rho^2 \mathbb{C}\_D}{\sqrt{\boldsymbol{k}}} \nabla \times \mathbf{b}, \\\\ \nabla (\boldsymbol{\varrho} \boldsymbol{v}) &= 0 \quad \boldsymbol{\omega} = \nabla \times \boldsymbol{v} \quad \mathbf{b} = \boldsymbol{v} \, |\boldsymbol{v}| \end{aligned} \tag{4}$$

where three more parameters are used to describe the properties of porous media: porosity *ϕ*, permeability *k* and drag factor *CD*. Porosity *ϕ* describes the volume ratio of pores to the total coil-filled sub-domain, *ϕ* = 1 − *Vcoil*/*Vsac*, where *Vcoil* is the accumulated volume of all coils, and *Vsac* is the volume of the aneurysm sac. The permeability *k* measures the fluid conductivity through porous media, *k* = *ϕ*3/ *cS*2, where *c* is the Kozeny coefficient related to the micro-shape of the porous media (for coils, the value of cylinders is chosen, *c* = 2), and *S* is the ratio of the surface area of all coils to the volume of porous region *Vsac*. The drag factor *CD* can be derived from the local Reynolds number. Note that when *ϕ* → 1 and *k* → ∞, these porous terms disappear, therefore, Equation 4 is identical to Equation 2, within the coil-free region. The computation of the extended porous terms takes little extra computational time (less than 1% of computing the other terms) when using the DEC method.

#### **5.2. Drag force**

In the existing simulations of aneurysm embolization, the interactive force between blood and coil was only studied for the blood from a global view, while the local reacting force on coils during the deployment process was ignored. In fact, the last term of Equation 4 is a description of the interactive force, but treated as an averaged quantity. When computing the reaction on the coil, we apply its local version, which is the drag force of flow over cylinder, since the coil is considered to consist of serially linked cylinder segments:

$$F\_D = \frac{1}{2} \mathcal{C}\_D \rho v\_\perp \left| v\_\perp \right| A \vert l\_\prime$$

where <sup>v</sup><sup>⊥</sup> is the velocity orthogonal to the coil, *<sup>A</sup>* is the cross-sectional area of the coil, *<sup>l</sup>* is the length of one short cylinder segment, and F*<sup>D</sup>* is the drag force applied on this segment. The velocity parallel to the coil is neglected, since it only produces shear force on the coil, which is insignificant compared to the drag force, and has little impact on the movement of the coil in the blood. Hence, the reacting force on the coil only depends on local fluid velocity.

#### **5.3. Fluid-structure coupling**

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

smoothly between the first two sub-domains. However, blood motion in all sub-domains is described uniformly by the Navier-Stokes equations (vorticity-based) of Brinkmann type:

<sup>Δ</sup><sup>ω</sup> <sup>−</sup> *μϕ*

∇(*ϕ*v) = 0 ω = ∇ × v b = v |v| ,

where three more parameters are used to describe the properties of porous media: porosity *ϕ*, permeability *k* and drag factor *CD*. Porosity *ϕ* describes the volume ratio of pores to the total coil-filled sub-domain, *ϕ* = 1 − *Vcoil*/*Vsac*, where *Vcoil* is the accumulated volume of all coils, and *Vsac* is the volume of the aneurysm sac. The permeability *k* measures the fluid conductivity through porous media, *k* = *ϕ*3/ *cS*2, where *c* is the Kozeny coefficient related to the micro-shape of the porous media (for coils, the value of cylinders is chosen, *c* = 2), and *S* is the ratio of the surface area of all coils to the volume of porous region *Vsac*. The drag factor *CD* can be derived from the local Reynolds number. Note that when *ϕ* → 1 and *k* → ∞, these porous terms disappear, therefore, Equation 4 is identical to Equation 2, within the coil-free region. The computation of the extended porous terms takes little extra computational time

In the existing simulations of aneurysm embolization, the interactive force between blood and coil was only studied for the blood from a global view, while the local reacting force on coils during the deployment process was ignored. In fact, the last term of Equation 4 is a description of the interactive force, but treated as an averaged quantity. When computing the reaction on the coil, we apply its local version, which is the drag force of flow over cylinder, since the coil

*CDρ*v<sup>⊥</sup> |v⊥| *<sup>A</sup>*|*l*,

*<sup>ρ</sup><sup>k</sup>* <sup>ω</sup> <sup>−</sup> *<sup>ϕ</sup>*<sup>2</sup>*CD* √*k*

∇ × b,

(4)

**Figure 9.** The shape of detachable coils after deployment.

*∂*ω

**5.2. Drag force**

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> <sup>L</sup>v<sup>ω</sup> <sup>=</sup> *<sup>μ</sup>*

*ρ*

(less than 1% of computing the other terms) when using the DEC method.

<sup>F</sup><sup>D</sup> <sup>=</sup> <sup>1</sup> 2

is considered to consist of serially linked cylinder segments:

For integrating the two models (coil and blood flow) in one single frame, we design a loosely-coupled strategy with the assumption that the simulation is performed over a series of identical cardiac cycles. Given the periodically time-varying boundary conditions at the inlet and outlet vessels around aneurysm, we solve the Navier-Stokes equations of Brinkmann type with a constant coil packing density by the DEC approach, and obtain the velocity at each tetrahedron center of the mesh at each time step. Then these velocity values are used to interpolate the velocity at the positions of coil segments and apply appropriate drag forces on the coil. The coil can provide real-time feedback inside the aneurysm at any time step during embolization. In this stage, we assume that a small segment of inserted coil does not change the blood flow. Until there is a significant increase of the coil density, the velocity of blood flow is recomputed at the new level of coil packing density. In this stage, we only care about the coil density, but not the coil shape. So the blood velocity field can be pre-computed and stored for real-time simulation of coil deployment. Although we use pre-computation of blood velocity, this process can be done in several minutes to simulate one cardiac cycle at five levels of coil density in our simulation. This is still essential for interactive planning.

The main benefit of the loosely-coupled approach is the relatively independency between the two systems, which allows different time resolution in two systems, independent real-time strategies, as well as pre-computation of the blood flow over one cardiac cycle. For the purpose of real-time refresh rate, we consider using relatively coarse mesh to reduce the size of the linear systems to be solved, and using large time steps to lessen the iterations necessary to simulate one second. As in other applications where real-time computation is sought, the objective is then to reach the best trade-off between accuracy and computational time.
