Numerical Simulation Model with Viscoelasticity of Arterial Wall

1. Introduction

1.1. Numerical analysis of intravascular flow

Numerical flow simulation is useful for understanding fluid phenomena such as blood flow or pulse wave propagation in the systemic arteries. For numerical analysis of intravascular flow, it is important to consider not only incompressible assumption and blood viscosity but also the viscoelasticity of the blood vessel wall; however, blood flow in vivo is complicated because of the unsteadiness of pulsatile flow and complex viscoelastic properties of the blood vessel wall. In order to conduct such numerical flow analysis in a viscoelastic blood vessel, it is effective to use the one-dimensional distributed parameter model, which can analyze flow along with the blood vessel axis. This distributed parameter model, pressure, flow volume and cross-section of the tube for every section element are defined and the time change is analyzed.

According to previous research, quantitative numerical simulation requires a model which take in both effects of unsteady viscous friction and viscoelasticity of the vessel wall in case flow unsteadiness is large (Reuderink et al., 1989). Conventional one-dimensional numerical simulation models can be classified into a linear distributed parameter model (Snyder et al., 1968; Avolio et al., 1980) and a nonlinear distributed parameter model (Anliker & Rockwall, 1971; Schaaf & Abbrecht, 1972; Porenta et al., 1986). The linear distributed parameter model has the feature that is easy to take in the influence of viscoelasticity and to conduct numerical analysis of the flow unsteadiness, since superposition of a periodic solution is possible; however, the influence of fluid nonlinearity cannot be disregarded. On the other hand, the conventional nonlinear distributed parameter model does not involve the effects of such flow unsteadiness and the viscoelastic behavior of the blood vessel wall concurrently with the difficulty of model construction. Hence, these models can be used only for qualitative discussion.

Consequently, in order to construct a numerical simulation model of intravascular flow for quantitative analysis with viscoelasticity of the blood vessel wall, it is necessary to use a nonlinear distributed parameter model to be able to include the viscoelasticity.

1.2. Clinical significance of blood vessel stiffness

Arterial stiffness indexes, such as PWV (pulse wave velocity: velocity of pressure wave propagation in circulatory systems), PP (pulse pressure: systolic blood pressure minus diastolic blood pressure), and AIx (augmentation index: index showing the effect of reflecting waves coming from peripheries) have received a lot of attention because they can indicate the risks for cardiovascular diseases (Oliver & Webb, 2003). These clinical indexes, which are obtained by analyzing blood pressure pulse waveforms in vivo, are used for predicting the prognosis of cardiovascular diseases and thus analyses of pulse waveform are clinically important. In order to understand the basis of these indexes, it is necessary to understand how changes in mechanical properties, such as blood vessel viscoelasticity, can affect pulse waveforms in the systemic arteries.

1.3. Research purpose and contents of this chapter

We thereby proposed a nonlinear one-dimensional numerical simulation model which can accurately calculate the viscous resistance of unsteady flow and the viscoelasticity of the tube wall. We have shown that the numerical model can describe well wave propagation in silicone tubes representing blood vessels (Kitawaki et al., 2003). The result showed that the viscoelasticity of the vessel wall plays an important role in the form of a pulsatile wave (Kitawaki & Shimizu, 2005); therefore, it is necessity to consider the viscoelastic effect accurately in the quantitative investigation of changes in pulsatile waves in vivo.

In section 2, we construct a one-dimensional numerical simulation model that takes into account unsteady viscosity and the generalized viscoelastic model from the theory of the fluid equation (Kitawaki et al., 2003; Kitawaki & Shimizu, 2006). In addition, the high-speed calculation method is described. When the wall of a blood vessel deforms finitely due to changes in the internal pressure, the wall’s physical properties, such as deformation compliance and viscoelasticity, change nonlinearly (Hayashi et al., 1980; Sato & Ohshima, 1985). In section 3, it was checked whether our one-dimensional model would be able to simulate the pulse wave propagation of small pressure waves in silicone tubes even when their deformation compliance and viscoelasticity changed independently, using appropriate values of the viscoelastic parameter of the silicone rubber tube for numerical simulation (Kitawaki & Shimizu, 2006). In section 4, the influence of viscoelasticity change on periodic pulsatile wave propagation was studied (Kitawaki & Shimizu, 2009). The purpose of the section was to investigate the effect of viscoelasticity change on periodic pulsatile wave quantitatively. For this purpose, we studied the pulse wave propagation of periodic small pressure waves using a silicone tube connected with terminal resistance, and obtained the waveform changes of pulsatile waves due to the change of mechanical properties, including the viscoelasticity of the tube. We then compared the experimental results with numerically calculated results using a one-dimensional numerical simulation model with terminal resistance treatment. Finally, we studied the effect of vessel wall viscoelasticity on the propagation of a periodic pulsatile wave by comparing numerical simulation results between the difference of viscoelastic models and viscoelastic parameters. The final section 5, describes the conclusion of this chapter.

2. Theory

2.1. Basic equations

When constructing the numerical model, we neglect the effects of bends of vessels. We also assume that the tube does not leak and that the flow is axisymmetric and incompressible. Under these conditions, the equations of continuity and momentum conservation of the one-dimensional model are given by (Olufsen, 1999),

$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=0$$E1

$$\frac{\partial Q}{\partial t}+\frac{\partial }{\partial x}\left(\frac{Q^2}{A}\right)+\frac{A}{\rho }\frac{\partial p}{\partial x}+F_t=0$$E2

where A is the cross-section of the tube, Q is the mean sectional flow volume, p is the mean pressure, t is time, x is distance along the vessel axis, ρ is the fluid density, and F_t is the viscous resistance.

By assuming that we are dealing with a Newtonian fluid, and an oscillating flow velocity distribution in a cylindrical tube, using the Womersley model, the viscous resistance F_t is given by (Zielke, 1968),

$$F_t=4\pi \nu \left\{2V+{\displaystyle {\int }_0^t\frac{\partial V}{\partial t}W(t-\tau )d\tau }\right\}$$E3

where V is the mean sectional velocity, ν is the kinematic viscosity, and W(t) is a weight function. In a rigid cylindrical tube,

$$W(t)=\underset{\omega \to \infty }{\lim }{\displaystyle {\int }_{-\pi }^{\pi }\frac{1}{F_J(\alpha i^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.})-2}{e}^{i\omega t}d\omega }$$E4

where$$F_J(z)=z{J}_0(z)/J_1(z)$$, and J₀, J₁ are 0th and 1st order Bessel functions of the complex number z. ω(=2πf) is the angular frequency of the flow oscillating at frequency f, $$\alpha =R\sqrt{\frac{\omega }{\nu }}$$is the Womersley number and R is the tube radius.

For a long wavelength, the flow velocity distribution in a distensible tube is similar to that in a rigid tube. Therefore, the weight function in a distensible tube can be approximated by Eq. (4), regardless of changes in the tube’s cross-section induced by the internal pressure. Note that the radius expressed as the Womersley number becomes a function of position and time, and consequently the weight function itself also becomes a function of position and time.

2.2. Deformation models of a wall

The tube law that describes the relationship between the tube cross-section and the internal pressure can be described by an elastic model and two viscoelastic models as below. In actuality, for a silicone tube, the relationship varies with internal pressure change, and shows nonlinearity. However, by assuming small local transformations, the following models can be used to approximate a linear model.

2.2.1. Elastic model

When the tube is perfectly elastic, the tube law of the elastic model can be expressed by,

$$\Delta p=\frac{1}{Cs}\Delta A$$E5

where Δp and ΔA are the differences in the pressure and cross-section relative to reference values, and $$Cs=\frac{\partial A}{\partial p}|\begin{array}{c}\\ A_0\end{array}$$ is tube deformation compliance at reference cross-section A₀.

2.2.2. Voigt model

By assuming that the viscoelasticity of the tube causes a phase lag between the applied pressure and resulting change in cross-section of the tube, from the Voigt model, the tube deformation law can be expressed by,

$$\Delta p=\frac{1}{Cs}\left\{\Delta A+{\tau }_V\frac{\partial A}{\partial t}\right\}$$E6

where τ_V is the relaxation time that accounts for the phase lag.

2.2.3. Generalized viscoelastic model

For complex viscoelasticity, as is the case for blood vessels, the generalized viscoelastic model can be applied (Westerhof and Noordergraaf, 1970). The following tube law of the generalized viscoelastic model can be derived, as shown in the next paragraph.

$$\Delta p=\frac{1}{Cs}\left\{\Delta A+{\displaystyle {\int }_0^t{\displaystyle \sum _{i=1}^n{f}_i\text{}\text{}{e}^{-(t-u)/{\tau }_i}\frac{dA(u)}{dt}du}}\right\}$$E7

The viscoelastic property of the tube wall is reflected in the second term of Eq. (7), which contains both the dynamic viscoelasticity parameter fi and the relaxation time parameter τ_i.

2.3. Generalized viscoelastic model

In order to obtain the tube law in one-dimensional flow analysis, the tube law of the generalized viscoelastic model is derived from the complex viscoelastic coefficient as follows. The complex viscoelastic coefficient of the generalized viscoelastic model E(s) can be expressed by (Westerhof and Noordergraaf, 1970),

$$E(s)=E_0\underset{i=1}{\overset{n}{\Pi }}(1+s\tau {\text{'}}_i)/\underset{i=1}{\overset{n}{\Pi }}(1+s{\tau }_i)$$E8

where E₀ is the static elastic coefficient, s is angular frequency (=iω), and τ_i and, τ'_i are relaxation time parameters which express the viscoelasticity of the tube wall. This equation can be transformed:

$E(s)=E_0\left(1+{\displaystyle \sum _{i=1}^n{f}_i}\frac{s}{s+1/{\tau }_i}\right)$E9

where

$$f_i=-\underset{k=1}{\overset{n}{\Pi }}({\tau }_i-\tau {\text{'}}_k)/\left\{{\tau }_i\underset{i\ne k}{\underset{k=1}{\overset{n}{\Pi }}}({\tau }_i-{\tau }_k)\right\}$$E10

Using this formula, the tube law in the frequency domain becomes

$$\Delta p(s)=\frac{h{E}_0}{2R{A}_0}\cdot \Delta A(s)\cdot \left\{1+{\displaystyle \sum _{i=1}^n{f}_i\frac{s}{s+1/{\tau }_i}}\right\}$$E11

where h is the tube wall thickness, Δp(s) and ΔA(s) are Laplace transformations of the pressure and cross-section from a reference value. The compliance of the tube deformation is$C_s=\frac{2R{A}_0}{h{E}_0}$.

By transforming Eq. (10) back to the time domain by inverse Laplace transformation, we obtain the tube law of the generalized viscoelastic model.

$$\Delta p=\frac{1}{Cs}\left\{\Delta A+{\displaystyle {\int }_0^t{\displaystyle \sum _{i=1}^n{f}_i{e}^{-(t-u)/{\tau }_i}\frac{dA(u)}{dt}du}}\right\}$$E12

2.4. High-speed calculation method

The following high-speed calculation method that was obtained for the viscous resistance term of rigid tube (Kagawa et al., 1983) was applied to the convolution integrals which appeared in the viscous resistance term in Eq. (3) and the viscoelastic terms in Eq. (7). By approximating the weight function in Eq. (4) by the sum of the exponential function,$$W(t)={\displaystyle \sum _{i=0}^k{m}_i{e}^{-n_i\tau }}$$; the viscous resistance in Eq. (3) can be expressed by recursive formulations as follows,

$$F_t(t)=4\pi \nu \left\{2V(t)+{\displaystyle \sum _{i=0}^k{y}_i(t)}\right\}$$E13

$$\{\begin{array}{ll}{y}_i(t)=0\hfill & (t=0)\hfill \\ \begin{array}{l}{y}_i(t+\Delta t)=e^{-n_i\Delta \tau }{y}_i(t)\hfill \\ \begin{array}{cc}& \end{array}+m_i{e}^{-n_i(\Delta \tau /2)}\{V(t+\Delta t)-V(t)\}\hfill \end{array}\hfill & (t>0)\hfill \end{array}$$E14

where τ=νt/R² and Δτ=νΔt/R² are normalized time and time step, respectively.

We need to determine the value of term number k in the Eq. (12) from the value of Δτ with consideration of approximation accuracy. In present experimental condition, the term number k was 10, because Δτ was calculated to be 1.2x10^-5.

The weight function of the convolution integral in the tube law of Eq. (7) can be expressed as a summation of the exponential function. Therefore, the tube law of the recursive formulations is as follows,

$$p(t)-p_0=\frac{1}{Cs}\left\{(A-A_0)+{\displaystyle \sum _{i=1}^n{z}_i}\right\}$$E15

$$\{\begin{array}{ll}{z}_i(t)=0\hfill & (t=0)\hfill \\ \begin{array}{l}{z}_i(t+\Delta t)=e^{-\Delta t/{\tau }_i}{z}_i(t)\hfill \\ \begin{array}{cc}& \end{array}+f_i{e}^{-\Delta t/2{\tau }_i}\{A(t+\Delta t)-A(t)\}\hfill \end{array}\hfill & (t>0)\hfill \end{array}$$E16

It is possible to determine the term number n in Eq. (14) from the viscoelastic characteristics.

3. Application to the fluid experiment of a computational model: In case of a difference in viscoelasticity

3.1. Method

3.1.1. Experimental model

Experimental apparatus

Figure 1 shows a schema of the experimental apparatus used in this study. The experiment tube consisted of two silicone tubes which was different form previous experiments (Kitawaki et al., 2003; Kitawaki & Shimizu 2005) with an inner diameter of 9 mm, a thickness of 0.5 mm, and a length of 1.45 m, connected by a rigid brass tube 5 cm in length and 9 mm in inner diameter. A piston pump was connected to one end of the tube via a 5-cm long brass tube, and a water tank with a valve was connected to the other end via another 5-cm long brass tube. The three brass tubes were fixed to a metal plate upon which the complete tube rested without longitudinal tension, allowing the silicone tubes to change shape freely. A pressure sensor (Nihon Koden, DX-100) was connected to the middle of each brass tube.

Figure 1.

Experimental apparatus

Figure 2(a) is a schema of the piston pump. The piston pump was driven by a computer-controlled stepping motor, and was capable of generating various waveforms with various flow volumes. The piston cylinder receives the backpressure from the water inside the tube. Therefore, displacement of the inner cylinder was measured by a laser displacement sensor (Keyence, LK-030).

Figure 2.

Piston pump

Experimental conditions

The tube, piston pump, and water tank were filled with water. Baseline internal pressures were set by adjusting the water head of the tank. The piston pump generated a single impulse. The flow volume at the inlet of the experiment tube was determined as shown in Fig. 2(b); when the impulse time t_i was set to 0.1 s, and the maximum flow volume Q_m was 5.0 ml/s. As seen in Fig. 2(c) which shows the Fourier transformation of the flow volume when a single impulse was generated, the highest frequency component of the flow volume was 30 Hz. Signals from the three pressure sensors and the displacement sensor were recorded for 10 seconds by a PC at a sampling rate of 1 kHz. The trials, that is, generation of a single impulse by the pump, were performed at increasingly higher baseline pressures of 2.5, 5.0, 7.5 and 10.0 kPa at 90-minute intervals. The trials were generated more than 70 minutes after changing the baseline pressures because 30 minutes were necessary for the viscoelasticity to return to a steady state. The baseline pressures were established by opening the valve, and the valve was closed just before the start of each trial.

3.1.2. Tube wall properties

Static tube law

Tube wall properties were determined for recent silicone tube. The static tube law and compliance of the silicone tube was obtained from the relationship between the volume of the piston pump and the internal pressure, as shown in Fig. 3, while performing one stroke of the piston pump over a period of about 30 minutes. The compliances, reference pressure, and reference cross-section of local tube deformation to the numerical calculation of each trial were determined from the collinear approximation of the pressure range of the experimental conditions. For example, a range of 2.5 - 4.3 kPa was used when the baseline pressure was 2.5 kPa, because the pressure pulse amplitude of a single impulse was about 1.8 kPa. These local compliances, assumed constant under the flow experimental conditions, are shown by squares in Fig. 3. These local compliances were used in all the wall viscoelastic models.

Figure 3.

Static relation between pressure and cross section of the silicone tube

Dynamic viscoelastic property of the tube

Dynamic viscoelasticity of a strip of the silicone tube was measured by a device (TA Instruments DMA2980). The measurement conditions were determined as initial strain of 0.18% and amplitude of 0.145%. These values correspond to a baseline pressure of 0.28 kPa and a wave amplitude of 1.8 kPa, respectively. The measured viscoelastic properties are in Fig. 4. This figure shows the dynamic modules of viscoelasticity normalized against the static modules in Fig. 4(a) and the loss tangent in Fig. 4(b), respectively. These results mean that, in the measurement range 0.1~30 Hz, the dynamic viscoelasticity property of the silicon tube, both the real parts and loss tangent, tends to increase gradually with increase in the frequency.

Figure 4.

Dynamic viscoelastic modules of silicone tube.(a) Ratio of dynamic and static modulus of viscoelasticity

The procedure used to decide the value of the viscoelastic tube parameter for the experimental condition change (Kitawaki & Shimizu, 2006) is described below. (1) First, the relaxation time parameter τ_i was determined so that it could cover the frequency range 0.1-30 Hz and it could increase between 0.01 and 200 Hz at regular intervals on a logarithmic scale. The term number n in Eq. (7) was set to 7 to keep the term number to a minimum and to show that the viscoelastic properties tended to increase smoothly. (2) The values of the dynamic viscoelasticity parameter fi were then determined to keep the difference between the experimental results and calculation results to a minimum from the low frequency term (i=1) to the high frequency term (i=7), using the different effect of each term on the numerical calculation results. (3) By changing the experimental conditions, the viscoelastic parameter was fitted by fixing the relaxation time parameter and changing only the dynamic viscoelasticity parameter.

Relaxation time τ_V of the Voigt model was determined as 2.5 ms from measurements of the delay time of the displacement of the outer diameter for the internal pressure change of the silicon tube.

3.1.3. Numerical calculation

The basic equations of the numerical calculation were digitized using a staggered grid system in space. For the calculation, Jameson-Baker’s 4th order 4 step method as a time differential and 4th order central differential with numerical friction as space was used (Jameson & Baker, 1983). Flow volume and cross-section of the next time step were obtained from an equation of continuity and momentum conservation, and then the pressure was calculated from the tube law as a function of time (Kitawaki et al., 2003). Convolution integrals appear in the viscous resistance term in Eq. (3) and in the viscoelastic term in Eq. (7). Since calculation of the convolution integrals requires a lot of computer memory to hold past velocity and cross-section values, and takes a lot of computational time, a high-speed calculation method for the viscous resistance term of a rigid tube (Kagawa et al., 1983) was applied, as shown in section 2.4.

The following conditions were used in the calculation.

1. Baseline pressure was determined as an initial pressure by the water head of the tank.

2. No flow in the initial state.

3. Flow volume calculated by displacement sensor was used as an input boundary condition.

4. No flow boundary condition was applied to the distal end.

Time step Δt and grid interval Δx were set at 0.5 ms and 0.05 m, respectively. The Courant number was 0.18~0.21 because the propagation velocity of the pressure wave was about 18~21 m/s, and the CFL condition (numerical stability condition) was satisfied. Actual calculation was performed on a workstation computer.

3.2. Results and discussions

3.2.1. Difference between calculated and experimental results

Experimental results of pressure propagation

The internal pressure waves measured at the 3 positions in the tube, and the flow volume are shown in Fig. 5, for a baseline pressure of 2.5 kPa. The flow volume was calculated from the measured displacement data of the piston pump by using LPF (FIR 25 Hz) and derivative filter.

Figure 5.

Time profiles of the pressure waves at three positions in the tube. The distance between input piston pump and measurement position are indicated on the left. The flow volume into the tube are given in the lower curve.

We can see that the pressure wave of the single impulse, generated by the movement of the piston pump, propagates towards the distal end, from where it is reflected. Upon returning to the proximal end, it is reflected again, back towards the distal end. The amplitude of the propagating pressure wave gradually attenuates, and the width of the pressure wave gradually increases because of the viscosity of the fluid and the viscoelasticity of the tube. The mean pressure rises for the fluid pushed by the piston pump. Both ends of the experimental tube are closed after the first pressure wave is generated by the piston pump. However, we can see an oscillatory wave when the pressure wave is returning back to the proximal end. This wave shows that the piston cylinder receives backpressure from the reflected pressure wave inside the tube.

Reynold's numbers and the wavelengths of the pressure waves were calculated from the measured value of the maximum velocity and pressure data from these experiments. Inspection of the Reynold's numbers (about 700) shows that laminar flow occurred under all experimental conditions. The wavelengths of the pressure waves (2.1 m) were sufficiently longer than the 4.5-mm radius of the tube, validating the long wavelength assumption and the assumptions of the Womersley model. Compared with the wavelength of the pressure wave, the length of the central rigid brass tube (5 cm) is sufficiently short, and because the rigid brass tube has a very small effect on the propagating pressure wave, it can be neglected in the calculation. The reproducibility of the pressure and flow volumes was good. Additionally, repeated measurements were virtually identical, so viscoelasticity changes such as memory effect did not happen.

Figure 6.

Calculated pressure waves (thin line) compared with the measured waves (bold line) at three location in the tube. Difference in viscoelastic models.

Comparison between calculated and experimental results

The numerical calculations were performed for a baseline pressure of 2.5 kPa and then compared with the experimental values of measurements of the pressure at the three positions along the experimental tube. The calculation results of combinations of the Womersley model with the elastic model, Voigt model and generalized viscoelastic model were compared, as shown in Fig. 6.

When the Womersley model was combined with all models, the initial pressure waves were close to the experimental value. However, when the Womersley model was combined with the elastic or Voigt models, the difference between the calculated and experimental results gradually increased because the calculated propagation velocities and attenuation level were underestimated. On the other hand, when combined with the generalized viscoelastic model together with the optimized viscoelastic parameter, the calculated results agreed well with the experimental results for all the experimental period. In the elastic model or Voigt model, when the tube deformation compliance or relaxation time τ_v was changed, the calculation result did not agree with the experimental result. These results show that the one-dimensional model using Womersley model combined with the generalized viscoelastic model are necessary in order for the numerically calculated result to agree with the experimental result.

3.2.2. Difference of viscoelasticity by the change of baseline pressure

Difference between pressure propagation experiments of baseline pressure difference

The flow volume of the inlet and the pressure waveform at the inlet position with changes in the baseline pressure are shown in Fig. 7. The flow volume in Fig. 7(a) are well controlled, and the movements of the piston are almost the same. As seen in Fig. 7(b), the pressure waveform changed differently when the baseline pressure was 2.5 kPa compared with other baseline pressures, because the propagation velocity was different. At baseline pressures of 5.0, 7.5 and10.0 kPa, the propagation is very similar during the first 0.8 seconds.

Figure 7.

Flow volume and pressure waves at beginning of the tube. Difference in initial pressures

Determination of viscoelastic properties

Because the experimental results for the baseline pressures of 5.0, 7.5 and 10.0 kPa are very similar, the optimized value of the viscoelastic parameter for baseline pressure of 5.0 kPa was used for the three experimental conditions, and comparisons of the calculated results using this viscoelastic parameter and the experimental results are shown in Fig. 8. As can be seen, the experimental and calculated results agreed well throughout the experimental period when the baseline pressure was 5.0 kPa. In contrast, at baseline pressures of 7.5 and 10.0 kPa, the agreement is good for the first 0.2 seconds, and gradually depreciates thereafter, because of the slight difference in their propagation velocities. Therefore, the experimental result for baseline pressure of 7.5 and 10.0 kPa cannot be simulated using the viscoelastic parameter optimized for baseline pressure of 5.0 kPa.

Figure 8.

Calculated pressure waves (thin line) compared with the measured waves (bold line) using viscoelastic parameter optimized for a baseline pressure of 5.0 kPa.

According to this result, the optimized value of the viscoelastic parameter for baseline pressure of 7.5 and 10.0 kPa were obtained respectively and calculations were conducted again. As shown in Fig. 9, the calculated and experimental results agree well. In the first 0.4 seconds, a clear difference between the two results can be seen, especially at higher baseline pressures. The differences could not be decreased by changing the viscoelastic parameters. Even with these differences, the agreement between the two results is good.

Figure 9.

Calculated pressure waves (thin line) compared with the measured waves (bold line) using optimized viscoelastic parameter for each experimental condition.

These results show that the one-dimensional model using the Womersley model combined with the generalized viscoelastic model can accurately simulate the effect of changes in the silicon tube’s viscoelasticity due to changes in the internal pressure when the viscoelastic parameter is appropriately determined. Additionally, the viscoelastic properties express the viscoelastic property change which depend the internal pressure of viscoelastic tube.

3.2.3. Relationship between viscoelasticity and static elasticity

The viscoelastic properties, both of the dynamic viscoelasticity parameters and the relaxation time parameters shown in Table 1 including the frequency f calculated from the relaxation time parameters, and the viscoelastic properties calculated from optimized values are plotted in Fig. 10. The squares in Fig. 10 indicate the measured values of the silicone tube viscoelasticity, and are close to the values obtained when the baseline pressure was 2.5 kPa, though a small error in the high frequency area of the loss tangent. This may be because the corresponding baseline pressure of the dynamic viscoelasticity measurement (0.28 kPa) is closer to the experimental baseline pressure of 2.5kPa than to the other experimental conditions.

From Fig. 10, it seems that the viscoelastic property change may be related to the baseline pressure. For example, when the baseline pressure increases from 2.5 to 5.0 kPa, the increase in the dynamic modules is steady at all frequencies. On the other hand, when the baseline pressure is 5.0 kPa or greater, the dynamic modules increase with increasing baseline pressure by increasing the frequency. Because the normalized dynamic modulus in the high frequency region changes with baseline pressure, the viscoelastic change can be said to be independent of the deformation compliance.

i	τi (s)	f(Hz)	fi
2.5kPa	5.0kPa	7.5kPa	10.0kPa
1	17.0	0.009	0.089	0.165	0.130	0.112
2	3.75	0.042	0.031	0.033	0.058	0.087
3	0.7	0.227	0.057	0.061	0.086	0.110
4	0.15	1.061	0.068	0.073	0.098	0.119
5	0.03	5.305	0.090	0.096	0.121	0.138
6	0.005	31.83	0.170	0.170	0.183	0.190
7	0.0008	198.9	0.220	0.220	0.223	0.224

Table 1.

Optimized viscoelastic parameters for each baseline pressure

Figure 10.

Change in dynamic viscoelastic module.

Accordingly, the one-dimensional numerical model, which takes into account unsteady viscosity and the generalized viscoelastic model, is good for simulating the propagation of small pressure waves in silicone tubes even when their deformation compliance and viscoelastic properties change independently.

3.3. Conclusion of this section

For the numerical analysis of the viscoelastic tubes, a nonlinear one-dimensional numerical model was investigated by including the unsteady viscous resistance and the effect of the tube wall viscoelasticity. By comparing the calculated results using these models with experimental results of a viscoelastic silicone tube, we can make the following conclusions.

1. The approximation error of the numerical simulation model is satisfactory small when the Womersley model is combined with the generalized viscoelastic model for the flow analysis of the viscoelastic tube even their deformation compliance and viscoelasticity change independently.

2. The relationship between the deformation compliance and viscoelasticity, which depends on the internal pressure, can be analyzed using this numerical model by appropriately choosing the value of the viscoelastic parameter.

4. Application to the fluid experiment of a computational model: In case of periodic pulsatile flow

4.1. Method

4.1.1. Experimental model

Experimental apparatus

Figure 11 shows a schema of the experimental apparatus used in this study in order to generate the periodic pressure pulse wave. A silicone tube with an inner diameter of 9.0 mm, 0.5 mm thick, and a 2.04 m long (L) was connected by three pressure sensors (Nihon Koden, DX-200), one each 2 cm from both ends of the tube and one in the center. These sensors were connected at intervals of 1.0 m. The experiment tube and each pressure sensor were connected by a narrow connecting tube made of silicone with an inner diameter of 3.0 mm, 0.5 mm thick, and 3 cm long. The frequency characteristics of each pressure sensor system, including these connecting tubes were sufficiently higher than those of the frequency component of the periodic pulsatile flow rate described later; hence, the sensor system can accurately measure the internal pressure of the experiment tube. A piston pump was connected to the proximal end of the experiment tube through the flow sensor (Nihon Koden, MFV-1100) and a tank was connected to the distal end of the tube through the terminal resistance. The experiment tube was placed on a metal plate without longitudinal tension, allowing the silicone tube to change shape freely. Figure 11(b) is a scheme of the piston pump connected to the proximal end of the experiment tube. The piston pump was driven by a computer-controlled stepping motor and was capable of generating various waveforms with various flow rates. Actual flow rate was measured at the proximal end of the tube by the flow sensor, as shown in Fig. 11.

Figure 11.

Experimental apparatus

The terminal resistance installed at the distal end of the tube simulated the peripheral arterioles in vivo. This terminal resistance, shown in Fig. 12, was filled with a bundle of around 200 thin stainless capillary tubes (26 G: with an outer diameter of 0.51 mm, inner diameter of 0.3 mm, and 50 mm long) in a rigid brass tube with an inner diameter of 8.0 mm.

Figure 12.

Terminal resistance

Experimental method

The experiment tube and piston pump were filled with water. Baseline internal pressure was set by adjusting the water head of the tank. Experimental trials were performed at increasingly higher baseline internal pressures of 5.6, 8.4 and 11.2 kPa (40, 60, 80 mmHg) respectively, because the deformation compliance and viscoelasticity of the experiment tube changed depending on the baseline internal pressure of the tube, as described later. The trials were generated at more than 60-minute intervals after changing the internal pressure, because 30-60 minutes were necessary for the viscoelasticity of the experiment tube to return to a steady relaxation state. Thus, the effect of a change in tube viscoelasticity on pulse wave propagation was examined.

Figure 13 shows the waveform of the flow rate (Q(t)) by the piston pump and the frequency component of the flow rate. In order to simulate periodic pulsation from the heart, we generated a pulsatile flow rate with 0.4-second ejection time (t_i) and about 2.8 ml/s maximum flow rate (Q_m) eight times in a 1.0-second period (t_p) at the proximal end of the experiment tube. In this case, the laminar flow condition was satisfied, as Reynold's number and Womersley’s number were 400 and 11.3, respectively. By the movement of the piston pump, the filling fluid from the piston pump flowed into the tank through the terminal resistance. The tank has a hole at the height of the baseline internal pressure head; thereby, the water head pressure always equaled the baseline internal pressure during each experimental trial. In each trial, signals from the three pressure sensors and the flow sensor were recorded for 10 seconds by a PC at a sampling rate of 1 kHz. The experiment was performed several times under the same conditions, and it was confirmed that reproducibility was high.

Figure 13.

Movement of piston pump

4.1.2. Numerical simulation

Mechanical properties of experiment tube

The static tube law and compliance of the experiment tube were obtained from the relationship between the increasing volume of the tube and the internal pressure while performing one slow stroke of the piston pump over a period of about 2.5 hours. This relationship between the static cross-section and the pressure, as well as the calculated tube deformation compliance are shown in Fig. 14. The static tube law was not exactly linear at any pressure. The measurement results in a range of pressures in each experimental trial were approximated to a linear model, because the nonlinear effects of tube law were not incorporated to the calculation in the present numerical simulation model. Tube deformation compliance at each baseline internal pressure is shown by squares in Fig. 14. At over 5 kPa, the tube deformation compliance increased with the increase in the internal pressure, showing countertrend to the tube law of blood vessels (Hayashi et al., 1980).

Figure 14.

Static relation between pressure and cross-section of the silicone tube

Viscoelasticity of a strip of the silicone tube was measured by a dynamic viscoelasticity measuring device. The measured viscoelastic properties are same as shown in section 3.1.2. This result indicates that both the dynamic modulus viscoelasticity ratio and loss tangent, tends to increase gradually with the increase in frequency. This trend is similar to the viscoelasticity property of blood vessels (Learoyd et al., 1966), and the dynamic modulus viscoelasticity ratio of the blood vessels is similar to or slightly more than that of the experiment tube. The procedure used to decide the values of the viscoelastic tube parameter are same as shown in section 3.1.2.

Terminal resistance

The viscous resistance of flow in the capillary tube of terminal resistance can be approximated by Poiseuille flow, as the terminal resistance occurred from rigid pipe flow and its inside diameter was sufficiently small. The ratio of the flow rate in the capillary tube to the pressure difference between both ends of the capillary tube (terminal resistance) seemed to be constant regardless of the flow rate; therefore, terminal resistance R_T (1.80 kPa/ml/s) was determined from the pressure difference and the flow rate obtained by moving the piston pump at a constant speed so as to generate a constant flow (about 1.1 ml/s).

Numerical simulation method

As initial conditions, baseline internal pressure and a cross-section which corresponded to the baseline internal pressure were used. It was assumed that there was no flow in the initial state inside the tube. Hence, the initial conditions can be shown as follows,

$$\{\begin{array}{c}p(x,0)=p_0\\ A(x,0)=A_0\\ Q(x,0)=0\end{array}$$E17

The flow rate measured by a flow sensor was used as the input boundary condition, and terminal resistance R_T was used as the output boundary condition. Hence, the boundary conditions can be shown as follows,

$$\{\begin{array}{l}Q(0,t)=Q(t)\\ Q(L,t)=\frac{p(L,t)}{R_T}\end{array}$$E18

We used the same numerical simulation methods as in a previous paper (Kitawaki & Shimizu, 2006), including calculation schemes, computational algorithms, and fast calculation methods. The basic equations of the numerical simulation were digitized using a staggered grid system in space. For the calculation, Jameson-Baker’s 4th order 4-step method as a time differential and 4th order central differential as space with numerical friction was used (Jameson & Baker, 1983). The flow rate and cross-section of the next time step were obtained from equations of continuity and momentum conservation, and then pressure was calculated from the tube law as a function of time. Convolution integrals appear in the viscous resistance term in Eq. (3) and in the viscoelastic term in Eq. (7). A high-speed calculation method was applied, since calculation of the convolution integrals requires significant computer memory to hold past velocity and cross-section values and requires much computational time.

Time step Δt and grid interval Δx were set at 0.5 ms and 0.04 m, respectively. The Courant number was 0.325 because the propagation velocity of the pressure wave was about 26 m/s at maximum. As a result, the CFL condition (numerical stability condition) was satisfied.

4.2. Results and discussions

4.2.1. Experimental results of pressure propagation

As an example of experimental results, the internal pressure measured at three points in the experiment tube, and the flow rate at the proximal end are shown in Fig. 15, at the baseline pressure of 5.6 kPa. The change in the whole time experiment is shown in Fig. 15(a), and an enlarged view of the last but one full-cycle pulsation is shown in Fig. 15(b). The pressure value shown here was set to be a differential pressure from the baseline internal pressure.

As shown in Fig. 15(a), there is a rise in the mean pressure or waveform change up to the first three pulses, but then it converges in the steady state. Also, there are oscillating pressure waves at the proximal and distal ends of the initial pulsation, due to wave reflections at both ends of the experiment tube; however, these oscillatory waves disappeared after the second pulsation, because this phenomenon is caused by the remaining reflected wave of the previous pulse.

As shown in Fig. 15(b) of the steady state of the pressure waveform, part of the waveform measured at the proximal end becomes almost flat in its upper portion for about 0.3 seconds, and then the pressure exponentially decreases. It is thought that this flat portion is generated by the progressive wave, which is produced by the movement of the piston pump, overlapping the reflected wave coming from the distal end of the tube. At the same time, the pressure waveform measured at the distal end has a pointed peak shape and a small reflected wave at the shoulder. The maximum pressure gradually increases from the proximal end of the tube towards the distal end, and the phenomenon is similar to the Peaking phenomenon in vivo.

4.2.2. Comparison between calculated results and experimental results

The calculated results in each tube law model were compared with experimental results in order to confirm the effect of the difference of the tube law model on the experimental results. Figure 16 shows examples of steady-state pulse waveforms at the baseline internal pressure of 5.6 kPa. Experimental results are shown by a bold black line, and calculated results are shown as a thin red line. In the Voigt model, calculated results in which the relaxation time parameter is multiplied by 4 (10.0ms) are also shown as a dotted line.

Figure 15.

Time profiles of pressure waves at three positions in the tube and flow volume into the tube.

In the calculated results using the elastic model and the Voigt model, the calculated maximum pressure is lower than that of the experimental results, and the calculated minimum pressure is higher than that of the experimental results; therefore, the difference between the calculated results and experimental results is large. Also, it is shown that the difference between the elastic model and the Voigt model is small, indicating that in the periodic change of waveforms, the Voigt model is not sufficiently effective. In addition, the calculated results did not agree with the experimental results in the Voigt model, because the effect is too small even when the value of the relaxation time parameter is varied as shown in Fig. 16(b).

Figure 16.

Calculated pressure waves (thin red line) compared with the measured waves (bold black line) at three locations in the tube

On the other hand, the calculated results using the generalized viscoelastic model almost completely agreed with the experimental results in the whole region. As mentioned above, it was proven that numerical simulation using the generalized viscoelastic model could express the experimental result accurately.

4.2.3. Difference in the viscoelastic parameter

In order to find out how many changes in the viscoelastic property of the tube were accompanied by baseline internal pressure change, and how this change influenced pulse wave propagation, the value of the viscoelastic parameter was determined under different conditions of baseline internal pressure. Initially, pressure waveforms at the baseline internal pressure of 8.4 kPa and 11.2 kPa were calculated using the value of the viscoelastic parameter at a baseline internal pressure of 5.6 kPa (hereinafter, this viscoelastic parameter is referred to as GVM(1)). The change in tube deformation compliance was incorporated into these calculations. As an example, a comparison of the calculated results and experimental results at the proximal and distal ends at a baseline internal pressure of 11.2 kPa is shown in Fig. 17. PP at the proximal end using GVM(1) was smaller than the experimental value, and it was proven that the viscoelastic effect of GVM(1) was not sufficient.

Figure 17.

Calculated pressure waves (dotted red line) compared with the measured waves (bold black line) at two locations in the tube

Based on this result, the value of the viscoelastic parameter was adjusted. We changed the value of the dynamic viscoelastic parameter in the low frequency region (i=1) so that the dynamic modulus elasticity ratio increased with the increase in baseline internal pressure in the experiment tube, and the value of the dynamic viscoelastic parameter in other frequency regions (i=2-7) was adjusted to a fixed ratio. The special evaluation function was not used to decide parameters but we took it into consideration to agree with the pulse waveform and the PWV. The viscoelastic parameter value decided in this way (hereinafter referred to as GVM(2)) is described in Table 2 and Fig. 18. The squares in Fig. 18 are the measured viscoelastic values of Fig. 4.

Agreement of the calculated results and experimental results using GVM(2), although not shown, was the same level as that in Fig. 16 (c).

Figure 18.

Change in viscoelastic properties of silicone tube

i	τi (s)	f(Hz)	fi
5.6kPa	8.4kPa	11.2kPa
1	17.0	0.009	0.095	0.169	0.227
2	3.75	0.042	0.033	0.034	0.036
3	0.7	0.227	0.061	0.063	0.066
4	0.15	1.061	0.073	0.076	0.079
5	0.03	5.305	0.096	0.100	0.104
6	0.005	31.83	0.170	0.177	0.184
7	0.0008	198.9	0.220	0.229	0.238

Table 2.

Optimized viscoelastic parameters for each baseline pressure

4.2.4. Effect of difference in viscoelasticity on pulse wave propagation

The effect of the difference in viscoelastic properties on pulse wave propagation was considered by analyzing the calculated results in which the viscoelastic model and the viscoelastic parameter were changed. PWV and PP were obtained from the calculated and experimental results. The relationship of these indexes with the baseline internal pressure is shown in Fig. 19. These indexes were obtained from mean values in four pulsations in the latter half in the steady state. Defining the maximum value of the second order differential waveform of the pulse waveform as an initial rise of the pulse wave, PWV was calculated by dividing the distance between calculation points (at the proximal and distal ends) by the time difference in the initial rise of these pulse waves. The error range of PWV is shown as an error in propagation velocity by a quantize error in calculating the propagation time. PP was obtained from the pulse wave at the proximal end.

Effect of dynamic elasticity modulus ratio

Figure 19(a) also shows the elastic tube theoretical velocity (Moens-Korteweg’s theoretical velocity) calculated from the tube inner diameter and tube deformation compliance. The theoretical velocity decreased by 5.8% as a result of the influence of the increase in tube deformation compliance, when the baseline internal pressure increased from 5.6 kPa to 11.2 kPa. On the other hand, PWV calculated from the experimental results was bigger than the elastic tube theoretical velocity, and the increasing ratio against the elastic theoretical velocity gradually increased with the rise in baseline internal pressure.

Figure 19.

PWV and PP change with internal pressure change

PWV using the elastic model approximately agreed with the theoretical elastic tube velocity, since the viscosity of the tube was not included. Using the Voigt model, the small effect by the viscoelasticity on the increase in PWV and PWV was slightly larger than the theoretical elastic tube velocity; however, the calculated and experimental results agreed well with the generalized viscoelastic model in GVM(2).

PP in the calculated results of the elastic model and the Voigt model was 90% of that in the experimental results. In the generalized viscoelastic model, the calculated results of GVM(2) agreed well with the experimental results.

This results show that, with the influence of viscoelastic change shown in Fig. 18, PWV increased by around 20-30% more than theoretical elastic tube velocity and PP increased by approximately 10% because increasing the elastic modulus by the effect of the low frequency component of viscoelasticity is more effective than the energy dissipation effect by the high frequency component of viscoelasticity. Furthermore, it can be said that the increase of viscoelasticity compensated for the influence caused by the change in tube deformation compliance.

Effect of the viscoelastic parameter

The effect of the viscoelastic parameter was examined from the difference in the calculated results between GVM(2) and GVM(1) for the baseline internal pressure of 11.2 kPa. The dynamic modulus elasticity ratio of GVM(1) was nearly 10.1% lower than that of GVM(2) within 1–10 Hz, as shown in Fig. 18. By this effect, PWV using GVM(1) was 6.6% lower than that using GVM(2), as shown in Fig. 19, and PP with GVM(1) was 6.4% lower than that with GVM(2). Considering the results with a baseline internal pressure of 8.4kPa, the difference of PWV and PP caused by the viscoelastic parameters was almost proportional to the change of baseline internal pressure; therefore, it might be necessary to determine the dynamic modulus elasticity ratio with accuracy below 4% in order to keep the accuracy of PWV and PP below 5%.

As described above, it was clarified that tube viscoelasticity played an important role in pulse wave propagation under periodic pulsatile conditions and had a significant effect on clinical arterial stiffness indexes, including PWV and PP.

4.3. Conclusion of this section

Using the one-dimensional numerical simulation model of a viscoelastic tube, the effect of tube viscoelasticity on periodic pulse wave propagation was examined by comparing the experimental results with the calculated results with a viscoelastic silicone tube. As a result, we can reach the following conclusions.

1. Numerical simulation using the generalized viscoelastic model can accurately express the experimental results with periodic pulsation of pressure waves.

2. When the dynamic modulus elasticity ratio increases, both PWV and PP increase, because increasing the elastic modulus is more effective than the energy dissipation effect by viscoelasticity change.

3. It is necessary to measure the dynamic modulus elasticity ratio with accuracy below 4% in order to estimate the values of PWV and PP with accuracy below 5%..

5. Conclusion

In this chapter, a nonlinear one-dimensional numerical model was constructed by including the effect of tube wall viscoelasticity and unsteady viscous resistance for numerical analysis of the viscoelastic tubes. This model can analyze the effect of viscoelasticity on intravascular flow that few previous papers have considered. Using a one-dimensional numerical simulation model of a viscoelastic tube, the effect of tube viscoelasticity on a pulsatile wave and periodic pulse wave propagation were examined by comparing the experimental results with the calculated results using a viscoelastic silicone tube.

Clinical blood vessel stiffness indexes vary with viscoelasticity change, because the elastic modulus is more effective than the energy dissipation effect by viscoelasticity change. This result showed that the viscoelasticity of the vessel wall plays an important role in the form of a pulsatile wave; therefore, it is important to consider the viscoelastic effect accurately in quantitative investigations using intravascular flow analysis.