Coupling of shear–circumferential stress pulses investigation through stress phase angle in FSI models of stenotic artery using experimental data

Abstract

Many cardiovascular diseases are closely associated with hemodynamic parameters. The main purpose of this study is mimicking a physiological blood flow in stenotic arteries to provide an understanding of hemodynamic parameters. An experimental setup was designed to produce original pulsatile flow and measure pressure pulse waves through a compliant tube. Moreover, a numerical model considering fluid–solid interaction was developed to investigate wall shear stress and circumferential stress waves, based on the results of the experiments. Results described elevated mean pressure by increasing stenosis severity especially at the critical obstacle of 50 %, which the pressure rose significantly and raised up by 10 mm Hg that may cause damage in endothelial cells. Increasing in stenosis severity led to: more negative wall shear stress and more oscillation of shear stress at the post-stenotic region and also more absolute value of angular phase difference between wall shear stress and circumferential stress waves at the stenotic throat. All of the aforementioned parameters determinant the endothelial cell pathology in predication of potential sites of progression of atherosclerotic plaques. Therefore, results can be applied in study of plaque growth and mechanisms of arterial remodeling in atherosclerosis.

Appendix: More details on the experimental setup

An overall view and schematic representation of the experimental setup for simulation of blood flow in elastic arteries is shown in Fig. 11. The setup consists of three major components: programmable pulsatile flow pump, elastic tube and the data acquisition and processing system.

Fig. 11

a Overall view of experimental setup, b comparison of flow waves as the output of the pump to the inlet of the test specimen

A custom-made pump was designed and manufactured to generate pulsatile flow in a wide range of arterial flow outputs. The device was comprised of mechanical and electrical units. The mechanical unit contained a servo-motor (MDFKS 056-23 190, Lenze, Germany), a planetary gearbox (MPRN 01, VOGEL, Germany), a ball screw (SFI2005, COMTOP, Taiwan) and a cylinder tank. The electronic unit contained a microcontroller (ATMega128, Atmel AVR^®, USA) for the control of the rotational pattern and speed of the servo-motor and consequent movement of the piston (Fig. 11). Due to the frequency response of the servo-motor (200 Hz) and the sampling rate of the microcontroller (1 m s), there were no limitations for producing flow pulses similar to those of human arteries. The inlet flow wave was produced through the computer interface based on previously published data of the brachial artery [24]. According to the considered flow wave, the maximum and minimum flow rates occurred at $\frac{t}{T} = 0.08$ and $\frac{t}{T} = 0.28$, respectively. Also, the mean flow rate was 4.92 (ml s^-1) [24].

To compare the flow wave produced by the pump to the flow wave at the inlet of the elastic tube, a rotameter with the range of 0–15 (ml s^-1) and a digital camera with accuracy of 80 (fps) were utilized. Results indicated negligible differences between the produced and measured flow, as demonstrated in Fig. 11b.

To study effects of arterial wall elasticity on blood flow parameters, an elastic tube with a defined stiffness modulus was used. For biological applications, medical grade silicon tubes have been used with differing dimensions and mechanical properties. In this study, an elastic tube (D-34209, B.Braun^®, Switzerland) with an inside diameter of 4.7 (mm) and wall thickness of 0.9 (mm) was used. This tube is comparable to an average human brachial artery with internal diameter and thickness of 4.5 (mm) and 0.83 (mm), respectively [40].

Mechanical properties of the tube were measured by a Universal Testing Machine (HCR 400-25, Zwick/Roell, Germany). The stress–strain relationship of the tube was similar to that of an arterial wall tissue with a hardening behavior, which is caused by elevation of pressure. Arterial walls within the human body experience pressure induced by circumferential strain pulse (diastole–systole) up to 10 % throughout the arterial tree [44]. For biological strain range, the elastic tube behaved with a linear stress–strain relationship. Hence, the stiffness modulus of the tube was considered as the slope of the curve within the aforementioned strain range. Results indicated elastic modulus of 463 (kPa) and Poisson’s ratio of 0.42 for the tube, which is similar to those of an average human brachial artery: 460 (kPa) and 0.42, respectively [40].

To produce stenosis, an external square metal housing and four bolts were used to adjust various area reductions (Fig. 2c). Due to the wall elasticity, significant thickness of the tube, and luminal pressure, it is assumed that the inner profile of the tube remains almost circular.

The data acquisition and processing system contained two pressure transducers (MLT0670, ADInstruments, Australia) that were used with the operational pressure range of 50−300 (mm Hg) (resolution of $\pm$ 0.1 (mm Hg)) and a processing unit connected to a computer. Each of the transducers was connected to a separate amplifier (ML117 BP Amp, ADInstruments™, Australia) in order to amplify the output signals before being processed. The amplifiers were connected to the main data processing unit (Powerlab/4SP, ADInstruments, Australia) by MLAC05 cables. The output pressure pulses were visualized using Chart for Windows™ v5.0.1 software in real time. The time lag between inlet and outlet pressure pulses was considered in evaluation of the pressure difference wave and is influential in calculation of shear and circumferential stresses. Such lag is due to pulse propagation within the distensible artery and depends on elastic modulus of the tube wall.

The experimental setup included elastic and resistant elements, placed before and after the elastic tube, respectively, as presented in Fig. 11. In order to resemble physiological data for the specific flow rate, we utilized the variable valve and the elastic element to regulate the pressure pulse and mean flow. In cardiovascular system, coupling of the heart to the aorta is of great importance. The highly distensible aortic root and ascending aorta play an essential role in this coupling, leading to a continuous aortic flow and proper function of the aortic tree. Such effect is described in the Windkessel model as the elastic element [44]. In order to adjust the amplitude of pressure pulses, an elastic element was used in between the pump and the elastic tube to control the pressure to make it comparable to diastolic–systolic pulse of a typical brachial artery. The mean value of the pressure pulse was tuned by a resistant element placed after the elastic tube. Such element represented the existing hemodynamic peripheral resistance of the distal circulation.

1.1 Governing equations

Governing equations of the fluid field: Continuity and momentum equations of Newtonian fluid with negligible body forces are described as:

$$ \left\{ {\begin{array}{*{20}l} {\uprho_{\text{f}} \partial_{t} \upnu_{\text{f}} + \uprho_{\text{f}} \nu_{f} \cdot \nabla \upnu_{\text{f}} - {\text{divT(}}\upnu_{\text{f}} ,\text{P}_{\text{f}} ) {\,=\,}C} \\ {{\text{div}}\upnu_{\text{f}} = 0} \\ \end{array} } \right. $$

(2)

which v_f (x,t) is the fluid velocity vector, P_f (x,t) is the fluid pressure, and ρ_f is the fluid density. The term ${\mathbf{T}}\left( {{\text{v}}_{\text{f}} ,{\text{P}}_{\text{f}} } \right) = - {\text{P}}_{\text{f}} {\mathbf{I}} + 2{{\upmu }}_{\text{f}} {\mathbf{D}}\left( {{\text{v}}_{\text{f}} } \right)$ describes the Cauchy stress tensor, in which μ_f describes the fluid viscosity and ${\mathbf{D}}\left( {{\text{v}}_{\text{f}} } \right) = \frac{1}{2}\left( {\nabla {\text{v}}_{\text{f}} + \left( {\nabla {\text{v}}_{\text{f}} } \right)^{\text{T}} } \right)$ is the deformation rate tensor.

Governing equations of the wall Large deformation theory was used to model arterial wall deformation. The continuity and momentum equations of elastodynamics for the wall are as follows:

$$ \left\{ {\begin{array}{*{20}c} {\uprho_{\text{s}} {\ddot{{\bf u}}} = \nabla .{\varvec{\upsigma}}_{\text{s}} +\uprho_{\text{s}} f} \\ {\frac{{\partial\uprho_{\text{s}} }}{{\partial {\text{t}}}} + \nabla .\uprho_{\text{s}} {\mathbf{u}}_{\text{s}} = 0} \\ \end{array} } \right. $$

(3)

where ${\mathbf{u}}_{\text{s}}$ is the wall displacement vector, $\sigma_{\text{s}}$ is the Cauchy stress tensor, ${\mathbf{f}}$ is the body force vector, and $\uprho_{\text{s}}$ is the wall density.

The Green–Lagrange strain was used due to large deformation of the wall. To do so, corresponding stress measure that relates stresses to the reference situation was adopted. The resulting symmetric stress tensor S, called the second Piola–Kirchhoff stress tensor, is a proper description for large deformation computations. To determine the stress, Cauchy stress has to be computed by:

$$ {\varvec{\upsigma}}_{\text{s}} = \frac{1}{\text{J}}{\mathbf{F}}.{\mathbf{S}}.{\mathbf{F}}^{\text{T}} $$

(4)

where ${\mathbf{F}} = \left( {\nabla_{0} {\vec{\text{x}}}} \right)^{\text{T}}$ is the deformation gradient tensor with respect to a reference configuration of the solid and J is the determinant of F. Thus, the Green–Lagrange strain tensor E could be defined as:

$$ {\mathbf{E}} = \frac{1}{2}\left( {{\mathbf{F}}^{\text{T}} {\mathbf{F}} - 1} \right) $$

(5)

Governing equations of the fluid–solid coupling: In order to obtain a complete FSI coupled system, two coupling conditions on the fluid–solid interface are required. The first condition describes the fluid–solid interface as a Dirichlet boundary for the fluid, i.e., the preset velocity values of the fluid must be equal to the structural nodal velocities, as is given by:

$$ {\mathbf{v}}_{\text{f}} = {\dot{\mathbf{u}}}_{\text{s}} \quad {\text{on T}}^{\text{c}} $$

(6)

where ${T}^{c}$ represents the interface between fluid and solid.

The second condition indicates that the fluid–solid interface is treated as a Neumann boundary for the solid, i.e., the solid boundary traction at the reference configuration, t _0s , is given by acting surface loads from the fluid boundary tractions, t _f , on interface at reference time as:

$$ {\text{t}}_{{0{\text{s}}}} = - \frac{{{\text{dT}}_{\text{ts}} }}{{{\text{dT}}_{{0{\text{s}}}} }}{\text{t}}_{\text{f}} \quad {\text{on T}}_{{0{\text{s}}}}^{\text{c}} $$

(7)

where dT is definition of differentiation due to load interchange on interface and subscripts 0 and t describe reference and changed configurations of fluid–solid interface, respectively. The negative sign refers to opposite directions of the normal vectors on the interface of fluid and structure with respect to a reference configuration. The formulation implies that the fluid boundary moves with the structure.

1.2 Numerical modeling

To further study effects of luminal stenosis on circumferential and wall shear stress waves together with the phase shift between them and the value of oscillation of shear stress, numerical models were analyzed, using fluid–structure interaction (FSI) method. The input data for models were based on experimental results. Models included degrees of severity of 20, 50 and 80 %. Models were developed and solved using FSI method by means of ANSYS CFX software as an axisymmetric model due to the symmetric geometry of an elastic tube and the applied stenosis within the experiment. Mechanical properties of the arterial wall (as the solid portion of model) and fluid field were allocated as described previously. The wall section was meshed with 8-node hexahedral mesh and the fluid field with 4-node tetrahedral mesh. Due to presence of stenosis, the mesh number was raised especially at stenosis site due to high stress gradient. A sensitivity analysis was performed for mesh number, and the optimized number was obtained in such a way that by elevation of mesh number, changes in resultant parameters were negligible. To verify grid independency, five grid sizes were applied for the model with 50 % stenosis. Table 1 shows the convergence of systolic average velocity at throat at systolic peak. A maximum difference of 1.6 % was observed when the density of 534 (elements mm^-3) was replaced by the density of 291 (elements mm^-3), and was assumed to be acceptable using the latter. For the solid domain due to the thin wall the density of 150 elements/mm³ was implemented.

Table 1 Average systolic velocity at stenosis throat for severity of 50 % for five grid sizes

Figure 2 shows part of a typical model of the stenotic artery with pre-stenotic and stenotic lengths of Z₁ = 30 (mm), Z₂ = 40 (mm) and total length of L = 110 (mm) for development of flow before and after luminal obstruction.

The applied flow wave [24] and the outlet pressure wave (P₂) measured and recorded in each set of experiments were used as two boundary conditions to the each respective model. The axial displacement of the arterial wall was restricted at inlet and outlet boundaries, as designed for the experiment. To reach a proper convergence, the model was pressurized to 80 (mm Hg) for 1 s and then the model was solved for five cardiac cycles with the defined period a cardiac cycle. To ensure convergence of results for all models, results of numerical models were calculated for the fifth cycle. This was achieved by sensitivity analysis of number of cycles based on minimal cycle to cycle change of results. Since resultant stresses (i.e., WSS and CS) for the fifth cycle were almost the same as the fourth cycle, we consider the fifth cycle for all reported stress results. The boundary conditions were applied far from the stenosis (5D from each end of the tube) to avoid the circumferential effects. Since for this region, the cross-sectional area is the same as the exit region, the additional pressure drops both in experiments and in numerical analysis were assumed to be negligible compared to that of the stenosis. The results were recorded from the fifth cycle at the mid-cross section of the arterial model to determine spatially averaged situation. Resultant parameters included inlet and outlet pressure pulses and pressure difference wave, wall shear and circumferential stress waves, oscillatory shear index and angular phase difference between WSS and CS in different scales of stenosis. The parameter OSI is defined as follows [5]:

$$ {\text{OSI}} = \frac{1}{2}\left( {1 - \frac{\left| \tau \right|}{\left| \tau \right|}} \right) $$

(8)

where $\left| \tau \right|$ is the magnitude of WSS and $\tau$ is the mean value of WSS.

The computational inlet pressure wave (P₁), compared to the P₁ wave measured and recorded in each set of experiment, was chosen as the criteria of validation of the model, as we chose appropriate P₂ and then solved the problem numerically to obtain P₁. Results showed a good agreement with both experimental and published data [24].

To analyze the phase difference between wall shear stress and circumferential stress waves a mathematical algorithm based on the discrete Fourier transform was carried. First we need to determine the waveforms of the diameter variation, D(t), and the wall shear stress, WSS(t), at a particular site of interest in the model. We decompose each waveform into a Fourier series with an amplitude and phase angle for each harmonic. The phase angle difference between the first harmonic of D(t) and WSS(t) is defined as the SPA. We can also define the SPA for the higher harmonics as we calculated in this study. Since a temporal variation of circumferential stretch is nearly in phase with that of pressure, we can approximate circumferential stretch by pressure variation [40]. The parameter SPA can be defined as follows:

$$ {\text{SPA }} = \varphi \left( {{\text{D}} - \tau } \right) = \varphi \left( {{\text{D}} - {\text{Q}}} \right) - \varphi \left( {\tau - {\text{Q}}} \right) \simeq \varphi \left( {{\text{P}} - {\text{Q}}} \right) - \varphi \left( {\tau - {\text{Q}}} \right) $$

(9)

Glossary

OSI

Oscillatory shear index. It is a parameter that shows us how much the wall shear stress changed in direction and magnitude (non-dimensional)

SPA

Stress phase angle. The phase angle difference between the first harmonic of pressure and wall shear stress (degree)