In Vitro Validation of Finite Element Analysis of Blood Flow in Deformable Models

Abstract

Introduction

The stress and strain in blood vessels, as well as hemodynamic parameters such as the three-dimensional blood flow and pressure fields, have direct effects on the initiation and development of cardiovascular diseases such as atherosclerosis and aneurysms.^{9, 10, 32} Knowledge of how in-vivo forces and tissue motions interact with implantable medical devices is also essential for understanding and predicting their behavior after implantation. For example, compliance mismatch between a prosthetic bypass graft and its adjacent native arteries has been hypothesized to lead to graft failure.²⁹ Medical imaging has been used to investigate vessel strain and blood flow hemodynamics, but with limited temporal and spatial resolutions, and often with discomfort to the patients as they are required to remain motionless for long periods during imaging sessions. Image-based computational fluid dynamics (CFD) methods, due to their minimal patient involvement and their ability to finely resolve time and space, have been a practical alternative to quantifying vessel strains and hemodynamic conditions for studies of disease mechanisms^{15, 31, 34} and the design and evaluation of medical devices.^{2, 4, 18} The ease of applying variations in geometry and flow conditions in the computational domain also motivates the use of CFD in the planning and prediction of surgical procedures.^{20, 30} Previous studies of cardiovascular CFD included the use of rigid wall models,²⁸ and dynamically deforming models.³⁷ Considering that vessel wall deformability often influences flow velocities and pressures, and that wave propagation phenomena can only be captured when considering wall deformation since blood behaves as an incompressible fluid, it is advantageous to include blood vessel deformability in numerical simulations whenever possible.

Much work remains to validate CFD methods against experimental data. Previous in-vitro validation studies have been performed only for the rigid case, likely due to the lack of realistic outflow boundary conditions, which are required to represent physical properties of downstream vasculature and to produce physiologic levels of pressure. For example, implementations of simple zero pressure boundary conditions in the physical setup where phantom outlets connect directly into a fluid reservoir^{1, 14, 16} would not be able to provide the pressures required to achieve realistic deformations in a compliant model.

In this study, we present results from two compliant phantoms under physiological flows and Windkessel boundary conditions for the validation of the numerical method incorporating wall deformability. The Windkessel boundary condition is a practical boundary condition prescription method in CFD simulations that can provide physiologically realistic impedances.^{11, 24, 33} We built a normal and a restricted physical model (flow phantom) comparable in size to the descending thoracic aorta, and constructed a physical analog of the Windkessel model to be attached to the outlets of each flow phantom to provide physiologically-realistic outflow impedances. A 1.5T MRI system was then used to acquire phase-contrast magnetic resonance imaging (PCMRI) flow velocity data in multiple 2D planes in the phantoms while they were under pulsatile physiologically-realistic flow and pressure conditions. The use of PCMRI in this study enables us to follow a similar protocol for future in-vivo validation. Next, we performed in-silico numerical CFD simulations incorporating a coupled momentum method for fluid-solid interaction (CMM-FSI) to include wall deformability,⁸ and prescribed Windkessel boundary conditions that directly corresponded to the physical setup. Flows, pressures, and velocity patterns measured in the in-vitro phantoms were then compared to those computed in the in-silico CFD simulations.

Methods

Physical Flow Phantom Construction & Characterization

We constructed two flow phantoms each containing a compliant vessel with an unpressurized diameter of 2cm, length of 25 cm, and thickness of 0.08 cm. The diameters and thicknesses of the vessels were selected such that physiological diameters^{13, 25} and strains^{6, 23} mimicking the descending thoracic aorta would be achieved under physiological operating pressures. One of the vessels was constructed to be a simple straight cylinder, and the other to be a straight cylinder containing a stenosis of diameter 0.88 cm at its center, which equates to an 84% area reduction relative to the mean operating diameter of the phantom.

To fabricate a compliant vessel, we used a multi-step dip-spin coating technique which entailed dipping an aluminum rod machined to the desired inner geometry of the vessel into a silicone mixture.^{3, 5} The important factors in obtaining the desired vessel wall thickness were the silicone mixture viscosity, dipping withdrawal speed, rod diameter, and number of dips. We set the silicone viscosity to 1500-2000 cp, dipped the aluminum rod vertically into the silicone mixture, and withdrew it at a controlled speed of 23.8 centimeters per minute. To obtain a uniform thickness circumferentially, we then set the rod on a horizontal rotating fixture for 30 minutes while the silicone dried. The entire process was repeated twice to obtain the desired thickness of 0.08 cm. Finally, the rod was set to cure in a heat convection oven at 100 °C for 4 hours.

We connected the inlet and outlet of each compliant vessel to a rigid section of flow conduit in order to allow for easy connection to the rest of the experimental flow setup. The bottom edge of the vessel was also glued along a ridge of width 6.8mm using a small amount of epoxy (5 Minute Epoxy, Devcon, MA) to mimic the in-vivo tethering of arteries to surrounding tissues such as the spine. The rigid inlet and outlet sections, together with the compliant vessel glued to the ridge, made up a “flow phantom.” For the stenotic phantom, a rigid ring was placed around the stenosis to make this region of the compliant vessel essentially “rigid”. This mimics the in-vivo material property of an arterial stenosis comprised of a stiff plaque.

Static pressurization tests were performed to characterize the phantom deformation under different levels of pressures. We pressurized the phantoms using a syringe while monitoring the internal pressure with a catheter pressure transducer (“Mikro-Tip” SPC-350, Millar Instruments, Houston, TX). The outer diameter of the compliant vessels at various pressurization levels was measured using a digital caliper (CD-6″ CS, Mitutoyo Corp., Kawasaki, Japan). Figure 1 shows the static pressurization characterization data for the two phantoms. There was an approximately linear relationship between the expansion of the compliant vessels and the increasing static pressures within the expected operating pressure range.²¹

Open in a separate window

Figure 1

Vessel Outer Diameter Versus Static Pressure for the a) Straight Phantom, and b) Stenotic Phantom

We characterized the viscoelastic properties of the silicone material using a dynamic mechanical analyzer (DMA Q-800, TA Instruments, New Castle, DE). Small rectangular pieces (20 × 5 × 0.7mm) of the silicone material were tested using a multi-frequency sweep controlled strain method. We applied a pre-load of 0.45N in order to deform the material to an approximately 12% static strain, we then applied an oscillating strain of 5% over a range of incremental frequencies from 0.10 – 6.0 Hz. These settings approximately correspond to the operating conditions of the phantoms in the flow experiments. When subjected to an oscillating strain at 1Hz, the storage moduli (which represent the elastic behavior of the material) in all of the six different samples tested were approximately 10% higher compared to when subjected to an oscillating strain at 0.1Hz.

Outlet Boundary Condition

A four-element Windkessel model consisting of an inductance (L), proximal resistance (R_p), capacitance (C), and distal resistance (R_d), was used at the outlet boundary of the phantom (Figure 2).¹⁷ Physically, the Windkessel module was designed such that physiologically realistic flows and pressures were achieved in the phantom, and that the specific values of the resistance and capacitance components remained reasonably constant over the operating range of each experiment.

Open in a separate window

Figure 2

In-vitro Flow Experiment Setup Diagram

We constructed the resistance module by placing a large number of thin-walled glass capillary tubes (Sutter Instrument, CA) in parallel with each other inside a plexiglass cylinder.¹⁷ Using Poiseuille's law and the equation for parallel resistances, the resistance value of the module is: Resistance=8 μ l/(Nπr⁴), where μ is the dynamic viscosity of the working fluid, l is the length of the capillary tubes, N is the total number of capillary tubes in parallel, and r is the inner radius of each individual capillary tube. We used 129 capillary tubes with inner diameters of 0.156 cm, and 125 tubes with inner diameters of 0.078 cm, for R_p and R_d, respectively. The theoretical resistance values corresponding to the specific viscosity of the working fluid in each of the two phantom experiments are listed in Table 1.

Table 1

Theoretical and Experimental Windkessel Component Values for the Straight and Stenotic Phantom Experiments.

	Straight Phantom Experiment		Stenotic Phantom Experiment
	Theoretical	Experimental	Theoretical	Experimental
L (Barye s2 cm-3)	7.0	7.0	7.0	7.0
Rp (Barye s cm-3)	250	240	240	260
C (cm3 Barye-1)	1.6 e-4	1.3 e-4	2.3 e-4	1.9 e-4
Rd (Barye s cm-3)	4100	4000	4000	4300

Open in a separate window

The capacitance of a fluid system is C=ΔV/ΔP where ΔV and ΔP are the changes in volume and pressure. In a closed system at constant temperature, an ideal gas exhibits the behavior PV=(P+ΔP)(V-ΔV), where P and V are the reference pressure and volume. The capacitance of a pocket of air is then: C = (V - ΔV)/P. For small changes in volume relative to the reference volume, a reasonably constant capacitance can be obtained with an air pocket. We constructed the capacitance module with a plexiglass box that can trap a precise amount of air to act as a capacitance in the system.¹⁷ The air volume used in the capacitor module was 210mL and 300mL at ground pressure (atmosphere pressure), for the straight and stenotic phantom experiment, respectively. The theoretical capacitance at the average operating pressure (relative to ground) of 100mmHg and 91mmHg for the straight and stenotic phantom experiments, respectively, is listed in Table 1.

The flow inductance results from the fluid momentum, and can be calculated from the geometry of the physical system. The flow inductance in a fluid conduit is: L = ρl/A, where l and A are the length and the cross-sectional area of the conduit, respectively. The theoretical inductance of the Windkessel module used in the phantom experiments is listed in Table 1.

In-vitro Experiment

We performed two in-vitro experiments, one with each flow phantom. For each experiment, the flow phantom and the outlet impedance module were placed in a flow system as shown in Figure 2. To produce the input flow to the phantom, we used a custom-built, MR-compatible, computer-controlled pulsatile pump to physically reproduce flow waveforms similar to the descending aortic flow measured in patients with aortic coarctations.²² Immediately upstream of the phantom, we placed one meter of straight, rigid tubing, a honeycomb flow straightener, and two pressure-stabilization grids, in order to provide sufficient entrance conditioning to generate a stable and fully-developed Womersley flow profile at the phantom inlet. The working fluid in the flow system was a 40% glycerol solution with a dynamic viscosity similar to that of blood, and contained 0.5% Gadolinium. The fluid viscosity for the straight and stenotic phantom experiments was measured to be 0.0461 poise and 0.0452 poise, respectively.

MR-compatible catheter pressure transducers (“Mikro-Tip” SPC-350, Millar Instruments, Houston, TX) were inserted through small ports on the sides of the phantom to capture pressure waveforms at various locations within the phantom (Figure 3), and also immediately downstream of the outlet impedance module. The signals from each catheter pressure transducer were sent into a pressure control unit (TCB-600, Millar Instruments, TX) which generates an electrical output of 0.5V per 100mmHg of pressure. An MR-compatible ultrasonic transit-time flow sensor was used to monitor the total input flow to the phantom. We placed the externally clamped flow probe (8PXL, Transonic Systems, NY) around a short section of Tygon tubing R3603 immediately upstream of the one-meter flow conditioning rigid tubing, and sent the signals from the probe into a flowmeter (TS410, Transonic Systems, NY) with its low pass filter setting at 160Hz. The data from the flow meter and the pressure control units was recorded at a sample rate of 96 samples per second using a data acquisition unit (USB-6259, National Instruments, Austin, TX) and a LabVIEW program (LabVIEW v.8, National Instruments, Austin, TX). For each data acquisition, approximately 50 cycles of flow and pressure data was averaged to obtain one representative cycle of flow and pressure waveforms. The flow and pressure waveforms were stable in between the cycles of each acquisition. We used the pressures measured downstream of the outlet impedance module as the ground reference, and subtracted it from all of the other pressure measurements to obtain the true pressure waveforms relative to the ground pressure.

Open in a separate window

Figure 3

Pressure and Flow Velocity Measurement Locations in a) the Straight Phantom, and b) the Stenotic Phantom. Green section is deformable. Grey section is rigid. Dimensions are in centimeters

We acquired through-plane flow velocity data at different locations within the phantoms (Figure 3) using a cardiac-gated 2D cine PCMRI sequence in a 1.5T GE MR scanner (Signa, GE Medical Systems, Waukesha, WI) and an 8-channel cardiac coil. The imaging parameters were: 256×192 acquisition matrix reconstructed to 256×256, 18×18 cm² field of view, 5mm slice thickness, TR=11∼14 ms, TE=5∼7 ms, 20 degree flip angle, and NEX=2. A velocity encoding gradient of 50 cm/s was used for all measurements, except for the L2, L3, and L4 measurements of the stenotic phantom, where the velocity encoding gradient was 100, 200, and 100 cm/s, respectively. The LabVIEW program which controlled the pulsatile pump produced a trigger signal that was converted by an electrocardiogram (ECG) simulator (Shelley Medical Imaging Technologies, London, Ontario, Canada) into an ECG signal used by the MRI system for gating, and 24 time points per cardiac cycle synchronized to this ECG signal were reconstructed. The temporal resolution of the velocity data was double the TR (∼26ms). We placed vitamin E capsules (Schiff Nutrition Group, Inc, Salt Lake City, UT) as well as saline bags around the flow region of each acquisition location to produce the reference signals of stationary fluids, which were then used for baseline eddy current correction with a linear correction algorithm in the analysis of the PCMRI data.

In-silico Simulation

We performed the numerical simulation of blood flow and pressure using a custom stabilized finite-element method to solve the incompressible Navier-Stokes equations.⁷ The deformability of the wall is incorporated by a CMM-FSI developed by Figueroa et al., which adopts a linearized kinematics formulation for the solid domain, and allows for a fixed fluid mesh and nonzero fluid velocities at the fluid-solid interface.⁸ The end result is that the effects of wall motion are embedded into the fluid equations simply as additional terms defined on the fluid-solid interface, leading to minimal increases in implementation complexity and computational efforts compared to rigid wall formulations.

Due to the fixed fluid mesh and linearized wall mechanics implementation of the CMM-FSI, we must use a mesh that most closely resembles the average geometry of the phantom during its operation. Calculations of the solid domain behavior require the definition of the elastic modulus of the vessel wall material. The data from the static tests of the compliant phantoms shown in Figure 1 can be used to find the radii of the compliant tubes at their respective average operating pressures, and the elastic modulus of the silicone material.

An analytical equation describing the expansion of a pressurized circular, cylindrical vessel made of an isotropic material and under small strain is:¹⁹

E=1.5ΔPRi2Ro(Ro2−Ri2)ΔR

(1)

where E is the elastic modulus of the material, R_i is the reference inner radius of the vessel, R_o is the reference outer radius of the vessel, ΔP is the change in pressure, and ΔR is the resulting change in radius.

Equation 1 describes a linear relationship between ΔR and ΔP. R_o is related to R_i by the vessel wall thickness, which we assume remains unchanged over small strains. Using the experimentally measured average operating pressure in the phantom, and prescribing the values of E and R_i, a plot of theoretical diameter versus pressure graph similar to Figure 1 can be generated. For each phantom, we modify the values of E and R_i until the theoretical plot coincides with the linear best fit of the static pressurization test results shown in Figure 1, and determine the operating radius of the phantom, as well as the elastic modulus of the silicone material. The operating radius was determined to be 1.13 cm at the operating pressure of 96 mmHg for the straight phantom, and 1.11 cm at the operating pressure of 89 mmHg for the stenotic phantom. The static elastic modulus of both phantoms was 9.1×10⁵ Pa. Using the result of the dynamic mechanical analysis previously discussed in the phantom characterization section, the effective elastic modulus of the silicone material at the fundamental operating frequency of the experiments (1 Hz) was then 1.0×10⁶ Pa.

At the outlet boundary, we prescribed a Windkessel model³³ (Figure 4) with the lumped-parameter component values determined from the experimentally measured flow and pressure at the phantom outlet. Mass conservation dictates that the average flow throughout different locations in the phantom is constant. For the analysis of Windkessel component values, we offset the outlet PCMRI flow measurement such that the average flow is equal to that measured by the flow probe, in order to exclude any effects of background correction variations. The measured average flow and pressure values at the phantom outlet were used to determine the total Windkessel resistance (R_p+R_d). The experimentally determined total Windkessel resistance in the straight and stenotic phantom experiment was 1.6% lower, and 7% higher, respectively, compared to theoretical. These ratios were used to scale the theoretical R_p and R_d to obtain the experimental values of the resistances. We then performed an analysis of the Windkessel model impedance:

Open in a separate window

Figure 4

Summary of Boundary Condition Prescriptions for the Numerical Simulations

Z(ω)=jωL+Rp+Rd1+jωCRd

(2)

to determine the value of the capacitance that will result in an impedance best reflecting the measured pressure and flow relationship at the phantom outlet. These experimentally determined Windkessel component values (as listed in Table 1) were then finally prescribed in the numerical simulations. For the stenotic phantom simulation, due to the presence of turbulence in the domain, an augmented Lagrangian method was used at the outlet to constrain the shapes of velocity profiles to prevent divergence.¹² This technique has been shown to have very little effect on the flow and pressure calculations in the numerical domain.¹²

We constructed the computational 3D solid models shown in Figure 3 from the physical construction details and the operating radii of each phantom. Each 3D solid model was discretized into an isotropic finite-element mesh with a maximum edge size of 0.1 cm using commercial mesh generation software (MeshSim, Simmetrix, Inc., NY). For the stenotic phantom mesh, we further refined a region of length 8cm distal to the stenosis using a maximum edge size of 0.03 cm, followed by another region of 4 cm downstream discretized using a maximum edge size of 0.06 cm. The straight phantom and stenotic phantom mesh contain approximately 1.5 million, and 3.5 million linear tetrahedral elements, respectively. Sections of the vessel wall boundary in the meshes were set to be rigid or deformable according to the physical construction of each phantom (Figure 4). We set the initial values of pressure and vessel wall distention in the mesh based on the average pressure in the physical experiment and the vessel wall properties. The flow waveform measured by the flow probe was then mapped to the inlet face of the computational domain using a Womersley velocity profile (Figure 4). A time step size of 0.42 milliseconds, which resulted in 2400 time steps per cardiac cycle, was used in the simulations. For the straight phantom simulation, we simulated 5 cardiac cycles and used the data from the last cycle where the pressures had stabilized in the final analysis. For the stenotic phantom, due to the presence of cycle-to-cycle variations in the velocity pattern, we simulated 14 cardiac cycles and used the ensemble-averaged data from the last 10 cycles in the final analysis. Since the PCMRI technique combines measurements acquired over multiple cycles into one cycle of velocity data, we ensemble averaged the stenotic phantom simulation results containing cycle-to-cycle variations to mimic the PCMRI data acquisition method.¹⁶

Results

Figure 5 compares the in-silico and in-vitro flow and pressure waveforms and normalized pulse amplitudes at various locations in the straight phantom. There is good agreement between the simulated and measured waveform shapes and amplitudes for both flow and pressure, throughout the various locations in the phantom. The average difference between the measured and simulated flows is 3.7 cc/s, which is 12% of the average flow (31 cc/s). The average difference between the measured and simulated pressures is 2.0 mmHg, which is 2.4% of the average pressure (85 mmHg). In both the simulation results and experimental measurements, we clearly observe progressive flow waveform damping, as well as progressive pressure waveform pulse amplitude increase, down the 25cm length of the deformable vessel. The most prominent pressure pulse amplitude increase occurs between the inlet and L1, which is the transition from the rigid section into the deformable section. The experimental measurements show that, of the approximate 10% increase in the normalized pressure pulse amplitude from the inlet to the outlet, roughly 7% occurred at the inlet and L1 transition. In addition, although not presented in Figure 5, we also observed a significant pressure waveform shape change between the inlet and L1. Lastly, both simulation and measurement showed an approximately 50% decrease in the normalized flow pulse amplitude between the inlet and the outlet of the phantom.

Open in a separate window

Figure 5

Straight Phantom Simulated Versus Measured Flow & Pressure a) Waveforms and b) Normalized Pulse Amplitudes, at Different Locations

Figure 6 shows the comparisons of flow and pressure waveforms and normalized pulse amplitudes at various locations in the stenotic phantom. The average difference between the measured and simulated flow waveforms is 3.4 cc/s, which is 13% of the average flow (26 cc/s). The average difference between the measured and simulated pressure waveforms is 1.0 mmHg, which is 1.2% of the average pressure (82 mmHg). We observe similar trends in pressure and flow behavior down the length of the stenotic phantom as those seen in the straight phantom. Between L2 and L3 of the stenotic phantom (across the stenosis), there is no visible difference in the flow waveforms, but there is a significant pressure drop. The drop in the peak pressure across the stenosis is 3.8 mmHg in the simulation, and 4.8 mmHg in the experimental measurement. The decrease in the normalized pressure pulse amplitude across the stenosis is 8.1% and 11.6% in simulation and measurement, respectively. A significant pressure waveform shape change occurred across the stenosis and is reflected in both measurement and simulation. As in the straight phantom case, the most prominent pressure pulse amplitude increase also occurs between the inlet and L1, which is the transition from the rigid section into the deformable section. The experimental measurements show that, of the 7% increase in pressure pulse amplitude between inlet and L2, about 6% occurred between the inlet and L1 transition. Both simulation and measurement show only an approximately 40% (compared to 50% in the straight phantom) decrease in the flow pulse amplitude between the inlet and the outlet of the stenotic phantom.

Open in a separate window

Figure 6

Stenotic Phantom Simulated Versus Measured Flow & Pressure a) Waveforms and b) Normalized Pulse Amplitudes, at Different Locations

Figure 7 compares the impedance modulus and phase at the phantom inlet and outlet between simulation and measurement. In both phantom experiments, there is agreement between the simulated and measured impedance across physiologically relevant frequencies. For the impedance phase, the experimental measurements exhibited more fluctuations at the higher frequency range compared to the simulated results. For frequencies in the 1∼3 Hz range, both the simulation and measurement show that the impedance modulus increases from the inlet to the outlet. For both phantoms, the general shapes and magnitudes of the impedance modulus and phase compare favorably with those measured in-vivo in previous studies.^{27, 35, 36}

Open in a separate window

Figure 7

Simulated Versus Measured Impedance Modulus and Phase at the Inlet and Outlet for the a) Straight, and b) Stenotic Phantom Experiment

We compare the simulated and measured through-plane velocity patterns at four different time points in the cardiac cycle: diastole, acceleration, peak systole, and deceleration. Figure 8 shows results for the L2 location in the straight phantom experiment. There is good qualitative agreement between the simulated and measured velocity pattern at all four time points. At acceleration and systole, there is forward flow of similar magnitudes and shapes across the slice, and a visible layer of decreased flow velocities near the vessel wall in both simulation and experiment. At diastole and deceleration, both simulation and experiment showed forward flow near the center, and a prominent region of backflow at the perimeter of the vessel. We generally observed a sustained Womersley velocity profile throughout the different locations in the straight phantom. Figure 9 shows flow velocity comparison results for the L2 and L3 locations in the stenotic phantom. In the L2 location, the pre-stenosis location in the stenotic phantom, we found nearly identical results as those presented for the straight phantom. In L3, the post-stenosis location, both simulation and measurement show a circularly-shaped jet of high forward velocities near the vessel center, and backward velocities (recirculation) around the vessel perimeter at the systole and deceleration time points. While the high velocity jet in the simulation contains slight irregularities in its shape, its size is comparable to that in the measured results.

Open in a separate window

Figure 8

Through-Plane Velocity Pattern Comparisons at the L2 Location for the Straight Phantom Experiment at Four Different Time Points: Diastole, Acceleration, Systole, and Deceleration

Open in a separate window

Figure 9

Through-Plane Velocity Pattern Comparisons at the a) L2, and b) L3 Location for the Stenotic Phantom Experiment at Four Different Time Points: Diastole, Acceleration, Systole, and Deceleration

Discussion

Figures 5 and ​and66 show that physiological flows and pressures^{21, 22} were achieved in the experiments. The damping of the flow waveform down the length of the phantom is the result of the flow being temporarily stored in the deformable tube, which essentially acts as a capacitance in the system. Since there is a rigid section at the center of the stenotic phantom, the total amount of deformable section is smaller compared to that in the straight phantom. With less deformable volume to absorb the flow pulse, the inlet waveform is better preserved and we indeed observed smaller flow waveform damping between the inlet and outlet in the stenotic phantom. The accurate numerical prediction of flow waveform shapes and magnitudes at different locations down the length of the phantoms indicates that the calculated fluid velocities at the vessel wall boundary accurately correspond to the physical vessel wall movement, faithfully capturing the compliant behavior of the vessel. The prediction of the decrease in the flow waveform peaks down the length of the phantom requires accurate calculations of the vessel wall expansion in response to the increase in pressure, where additional fluid is stored in the vessel during the systolic period. On the other hand, the prediction of the gradual increase in flow waveform minimums requires accurate calculations of the elastic behavior of the vessel wall, which releases the stored fluid during the diastolic period. The flow waveform damping behavior observed in the simulations and experiments is generally consistent with the actual blood flow behavior in-vivo, where the pulsatility resulting from the pumping heart is damped out throughout the vasculature, and eventually transformed into steady flow in the capillaries.

Since a significant portion of the cardiac cycle was diastole where the flow was low or retrograde, the average flow rates were relatively low compared to systolic flow rates. The average difference of approximately 3.6 cc/s between the measured and simulated flow was 12∼13% relative to the average flow rates, but was only under 5% relative to the systolic flow rates, which were between 80∼140 cc/s.

The progressive increase of pressure pulse amplitude down the length of the phantom is also consistent with the in-vivo observation where the pulse pressure progressively increases from the brachial artery downstream towards the radial artery.²⁶ It has been generally believed that such phenomenon is attributed to the increased stiffness of the downstream arteries. Our experimental and simulation results suggest that wave propagation and reflection alone could contribute to a pressure pulse increase under the condition of constant vessel stiffness.

Across the rigid and deformable junction where there is a mismatch in characteristic impedances, the change in the pressure waveform shape and the prominent increase in the pressure pulse amplitude could be attributed to wave reflections. Across a stenosis, pressure waveform changes also occur due to energy losses in the post-stenosis turbulent flow region. The numerical simulation accurately captured both the wave reflections between the rigid and deformable sections, and the energy loss across a stenosis, accurately predicting the changes in the pressure waveform at different locations within the straight and stenotic phantoms.

The impedance modulus increase between the phantom inlet and outlet during the 1-3Hz frequency range reflects the capacitive effect of the deformable tube. Pulsatile flow enters the inlet with relative ease due to the compliance in the deformable vessel downstream, resulting in a lower impedance modulus. At the outlet of the phantom, there is no deformable region downstream to manifest the effect of the lowered resistance to pulsatile flow. This phenomenon is only prominent in the lower frequency range partly because of the physical characteristics of the deformable tube dictating its response to dynamic strain, and partly because of the low modulus values in the higher frequency region making any differences difficult to observe. The increased prominence of outlet impedance phase oscillations at the higher frequency region in the experimental measurements could be due to the small high frequency component in the flow and pressure waveforms, making the noise in the measurements relatively high.

The favorable comparison in velocity patterns between simulation and measurement for the straight phantom is consistent with our expectation due to the trivial geometry of the phantom. We also expect the pre-stenosis location in the stenotic phantom to show similar results since complex flow originates from the stenosis and propagates to the regions downstream. At the L3 location in the stenotic phantom, which is immediately downstream of the stenosis where complex and recirculating flow occurs, the simulation showed a smooth contour for the flow pattern right up to the time frame immediately prior to peak systole, after which point the flow begins to decelerate and diverge, resulting in slight irregularities in the shapes of the flow patterns in the simulation results. We found that the irregularities in the flow pattern shapes were correlated with mesh resolution, where a simulation computed on a finer mesh resulted in fewer irregularities. In a previous simulation we performed using a mesh without local mesh refinements in the stenosis and its downstream regions (resulting in a mesh containing approximately 1.5 million elements), we observed similar irregularities that were much more pronounced. In regions containing complex and diverging flow, it is thus important to define a desired balance between flow pattern prediction accuracy and computational cost.

The vessel wall motion in the numerical simulation is sensitive to the prescribed thickness and elastic modulus of the vessel wall. Both the wall motion and the prescribed vessel geometry affect the volumetric flow and the pressure changes down the length of the vessel. Prescription of higher elastic modulus or smaller vessel diameter would result in diminished wall motion and smaller flow waveform damping, and vice versa. The method we developed to determine the relevant geometry and vessel wall properties requires direct manipulations and observations of the vessel which are only possible in-vitro. To apply the numerical simulation in an in-vivo setting, additional methods would need to be developed to determine the equivalent values of vessel parameters for the numerical model.

In conclusion, in this study we have produced a set of in-vitro, high-quality experimental data that can be used to compare against CFD results of flow and pressure within a compliant vessel under physiological conditions. The deformable CFD simulation utilizing the CMM-FSI and a fixed fluid mesh was capable of capturing realistic vascular flow and pressure behaviors. There were good predictions of flow and pressure waveforms down the length of a straight and a stenotic deformable phantom, indicating that the numerical simulation captured both the vessel wall motions and wave reflections accurately. Due to the good comparisons in pressure and flow, the impedance comparisons were also favorable. The simulated and measured flow and pressure results were similar to those previously measured in-vivo. The numerical simulation was able to track velocity patterns very well in regions with simple flow. In regions containing more complex and diverging flow, a finer mesh resolution was required for the simulation to capture the velocity patterns faithfully. The results presented in this paper show promising potential for the numerical technique to make accurate predictions of vascular tissue motion, and blood flow and pressure under the influence of blood vessel compliance. This study provides the cornerstone for further deformable validation studies involving more complex geometries, and in-vivo validation studies that could ultimately support the use of CFD into clinical medicine.

Acknowledgments

References