## Abstract

Cardiovascular diseases, mostly related to atherosclerosis, are the major cause of death in industrial countries. It is observed that blood flow dynamics play an important role in the aetiology of atherosclerosis. Especially, the blood velocity distribution is an important indicator for predisposition regions. Today magnetic resonance imaging (MRI) delivers, in addition to the morphology of the cardiovascular system, blood flow patterns. However, the spatial resolution of the data is slightly less than 1 mm and owing to severe restrictions in magnetic field gradient switching frequencies and intensities, this limit will be very hard to overcome. In this paper, constrained fluid dynamics is applied within the smoothed particle hydrodynamics formalism to enhance the MRI flow data. On the one hand, constraints based on the known volumetric flow rate are applied. They prove the plausibility of the order of magnitude of the measurements. On the other hand, the higher resolution of the simulation allows one to determine in detail the flow field between the coarse data points and thus to improve their spatial resolution.

## 1. Introduction

Cardiovascular diseases related to atherosclerosis are the most common cause of death in industrial countries. Observations are made that haemodynamics play an important role in the aetiology of atherosclerosis [1,2]. Especially, the wall shear stress and consequently the velocities of the blood flow are important indicators, because they turn out to be smaller in the predisposition regions [3,4].

In this context, blood velocity distributions in the cardiovascular system can help detect atherosclerosis at very early stages. Phase-contrast magnetic resonance imaging (MRI) is a flow-sensitive imaging technique that allows not only the non-invasive depiction of human body morphology but also the acquisition of time-resolved blood flow velocities in three directions [5,6]. Nevertheless, the quality of the flow-sensitive data remains challenging. The revealed datasets do suffer from a limited spatial resolution of slightly below 1 mm due, among others, to severe restrictions in magnetic field gradient switching frequencies and intensities. The rapid switching of the magnetic field gradients can result in nerve stimulation and limitations in gradients’ intensities are related to instrumentation.

As we will show, numerical modelling and, especially, the implementation of constraints provide an approach to enhance the measurements. A constraint of the form , where denotes the particles’ velocities, is an additional condition that has to be solved together with the equations defining the physical system. In this paper, a new method for the implementation of velocity constraints within the smoothed particle hydrodynamics (SPH) formalism is presented.

The new method is first applied to the validation of MRI flow data measured in a physical phantom. The constraint is based on a volumetric flow rate known from the experiment. Furthermore, local velocity constraints are introduced. In this case, the velocities measured in a pathological model are imposed on the corresponding simulation regions. Thus, a more detailed observation of the velocity distribution with improved spatial resolution is possible.

## 2. Magnetic resonance imaging flow data

In this work, experimental data provided by the group of Dr Michael Markel from the University Clinic in Freiburg, Germany were used [7]. The objective of the experiments is to model parts of the human aorta with simple geometries. For example, the descending thoracic aorta is modelled with a straight rigid polyvinyl chloride tube. The tube is connected with two conical adapters and two hoses to a pump in order to produce a constant flow. The resulting flow is measured using the phase-contrast MRI method. The fluid used is a glycerin–water 30 : 70 vol% mixture, which is known to have a viscosity behaviour similar to that of blood.

The resulting flow dataset has a resolution of 1.0417×1.0417×1 mm^{3}, resulting in 35 voxels. The extracted velocity profile from the middle part of the tube is shown in figure 1*a*. The unfitted data do not show a perfect parabolic profile expected for the laminar flow in a straight tube. Deviations stem from noise, measurement inaccuracies and deviations from ideal laminar flow behaviour.

In the same figure, the resulting profile from an SPH simulation with constrained volume flow rate is shown. The simulation set-up is explained in §6. It can be seen that the constrained simulation delivers the expected analytical profile. Indeed, for the tube we might deduce this from the analytical solution but for more complex geometries we require the SPH simulations.

## 3. Smoothed particle hydrodynamics

### (a) Equations of motion

The fundamentals of SPH are extensively presented in earlier studies [8–11]. For the simulation of viscous incompressible flow, the equations of motion are the Navier–Stokes equations, describing conservation of momentum and the continuity equation describing mass conservation. In their Lagrangian form, the equations read
3.1and
3.2where **v** denotes the velocity of the fluid, *P* the pressure, *ρ* the density and *η* the dynamic viscosity of the fluid.

In the following, the indices indicate the respective SPH particles and *f*_{i} is the value of the function *f* to be taken at the position of particle *i*. The ∇*P*/*ρ* term is discretized as proposed in the study of Monaghan [11], thus producing an anti-symmetric expression that conserves linear and angular momentum exactly:
3.3The most common way to discretize the viscous term is [8]
3.4Within the quasi-incompressible approach [12], the continuity equation allows for density gradients and is given by
3.5Discretizing as proposed in the study by Monaghan [8], one obtains
3.6where **v**_{ij}=**v**_{i}−**v**_{j}, **r**_{ij}=**r**_{i}−**r**_{j}. Furthermore, the scalar function *w*_{ij}=*w*(|**r**_{ij}|) was introduced. It is defined by ∇*W*_{ij}=−**r**_{ij}*w*_{ij}.

An artificial equation relating the local density *ρ* to the pressure *P* is the Tait equation. It limits the density fluctuations sufficiently and is defined by
3.7where *γ*=7, *ρ*_{0} is a reference density and *P*_{0}=*c*^{2}*ρ*_{0}/*γ*. For the reference density, the actual density of the material is used.

*P*_{0} is determined by the speed of sound *c*. The speed of sound *c* is estimated as proposed in Morris *et al*. [13] to keep the density fluctuation below 3 per cent.

### (b) Correction of boundary deficiencies

SPH suffers from kernel sum deficiencies at boundaries. The summation relies on a complete neighbourhood of the surrounding particle *j* for the particle *i*, for which a property should be calculated. This is only the case for the particles in the bulk. In this work, a solution presented earlier [14,15] is applied. A corrective version of the kernel is obtained involving the corrective matrix **B**, which is defined as
3.8where *f*_{ij}=*f*_{i}−*f*_{j}.

Thus, the gradient of a scalar field *f* can be computed using
3.9

### (c) Imposing no-slip conditions

No-slip boundary conditions are realized by regarding the viscous interactions **F**_{η}(**v**_{i},**v**_{j}) between fluid particles and the so-called wall particles. Detailed description can be found in Morris *et al.* [13].

## 4. Velocity constraints

In the current paper, a new method for the implementation of velocity constraints is presented. It does not include iterative procedures, but the linear condition can be solved exactly using a predictor–corrector time-integration scheme. The scheme reads
4.1
4.2
4.3
and
4.4
where Δ*t* is the time step, **p** the momentum, the predicted value for the momenta at time *t*+Δ*t* and **F** the force on the particle *i*. The system is first advanced over a single time step. Ignoring the constraints, the resulting momenta from the prediction–correction algorithm are . Then they must be adjusted to **p**_{i}(*t*+Δ*t*) to satisfy the velocity constraint
4.5
where **v**_{k} is the mean velocity we want to impose, e.g. the volumetric flow rate which is measured by an experiment.

The correction is calculated using the Lagrange multiplier formulation [16]:
4.6
From equations (4.5) and (4.6), it follows that the velocity correction **v**_{corr} for particle mass *m*=1 applied at the end of the integration must be
4.7
In the next time step, the corrected velocities are used for the computation of the new positions of the particles.

## 5. Validation

The validation of the SPH formalism with constraints is based on a Poiseuille flow between two parallel plates. In this case, the analytical solution for the flow in the *z* direction and walls in the *y* direction is a parabolic profile given by [16]
5.1
where *F*_{z} is a homogeneous body force and *D* the distance between the plates. The walls are sitting at *x*=0 and *x*=*D*.

It is expected that the force resulting from imposing the velocity constraint (4.5) with *v*_{k}=*v*_{z,mean} would be equal to the body force *F*_{z} producing the same velocity profile. To stay close to the experimental data in this simulation, the same fluid viscosity *η*=3×10^{−3} Pa s was used. During the simulation with constraints, first the correction velocity *v*_{z,corr} (4.7) is measured. The resulting force is calculated by
5.2

In figure 1*b*, the homogeneous correction force *F*_{z,corr} on the particles is plotted over the simulation time. The initial velocity conditions are a parabolic velocity profile. Therefore, the correction in the first time step is zero. The curve increases up to an overshoot. This overshoot is caused by numerical deviations between the instantaneous mean velocity, , of the particles and the imposed velocity *v*_{k}. The uncertainty vanishes within a few time steps. The saturation value of the force is 4.25 N. Compared with the analytical value *F*_{z} of 4.43 N, the deviation is roughly 4 per cent. This is of the order of the error we allow for the density fluctuations. Hence, these values are approximately equal. Thus, we have verified that the body force and the formalism with constraints do yield the same result.

## 6. Simulations with constraints

The first application of the SPH formalism with constraints is to prove the plausibility of the measured MRI data regarding their order of magnitude.

From the experiment presented in §2, the volumetric flow rate in the physical phantom (tube) is known. It corresponds to a mean velocity of 14.62 cm s^{−1}. An SPH simulation with the constraint *σ*_{k} (4.5) and **v**_{k}=14.62 cm s^{−1} was performed in a periodic box with walls in one direction. The simulated viscosity corresponds to that of the experiment, *η*=3×10^{−3} Pa s, and the box length to the tube diameter. Thus, a Reynolds number of 1600 is simulated. The distance between the walls is discretized by 70 particles, thus producing twice as many measurement points as in the raw MRI data and thus increasing the resolution from 1.0417 to 0.52 mm. The resulting profile is shown in figure 1*a*. It can be seen that it matches with the data fit and thus verifies the order of magnitude of the measured MRI data. This verification based on the SPH formalism with constraints can be applied to arbitrarily complex geometries of blood vessel models as soon as the measured volumetric flow rate is additionally known.

Another application of the method for the implementation of constraints is the fast modelling of different pathologies in the cardiovascular system. The two datasets in figure 2*a*,*c* correspond to one-side plaque and a stenosis. The regions pointed with the arrows are the deposition regions. The measured velocity in these voxels is **v**_{0}=0.

These anatomical situations can be approximated by a periodic box with walls in one direction. The fluid viscosity is again *η*=3×10^{−3} Pa s corresponding to the fluid used for blood modelling in the MRI measurements. Twenty particles are created between the walls, thus giving twice as many points as the 10 voxels in the dataset. The pathology regions are constrained according to the measurements. The aim is to observe the velocity field resulting from the simulation, which shows a resolution higher than that of the measured data. These are shown in figure 2*b*,*d*, where the particles’ velocities are represented vectorially. In both figures, it can be seen that the particles in the middle show the highest velocities. Thinking about a flow going into a narrower tube, this is a natural result. Behind the constrained region a vortex is observed. This indicates a predisposition region for further accumulation of plaque.

## 7. Summary and discussion

This work presents a method for the implementation of velocity constraints within particle methods. With the introduced method, an SPH enhancement of MRI flow data can be achieved. It avoids complex iterative procedures and can be easily implemented in the time-integration algorithm.

The constrained SPH method is used for the validation of MRI flow data. Based on a known volume, flow rate constrained simulations can verify the order of magnitude of the measurements. In this work, we concentrated on MRI flow data measured in a simple physical phantom, but the method can be extended to more complex geometries. Furthermore, the new method for the implementation of constraints can be used for the fast modelling of cardiovascular disorders, like stenosis. The finer resolution of the simulation allows a deeper insight into the velocity distribution and can reveal predisposition regions for plaque accumulation.

Several directions are worth exploring for future investigations. A complete aorta model can be simulated and the SPH formalism used can be expanded to model blood flow phenomena such as non-Newtonian rheology, moving arteries and complex fluid–wall interactions. However, only limited knowledge of the physics of these phenomena is available. Finally, the particles can be used to model each blood component by a discrete particle or, alternatively, a continuous concentration can be assigned to each fluid particle as an additional degree of freedom [17,18].

## Acknowledgements

The authors are grateful to the University Clinic in Freiburg, Germany for providing measurements data and the German Science Foundation DFG for financial support in the framework of SFB499. D.K. and J.G.K. acknowledge funding by the University of Freiburg through the German excellence initiative.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society