Verification of fluid-structure interaction modelling for wave propagation in fluid-filled elastic tubes

Abstract

This paper presents a verification study of wave propagation in fluid-filled elastic tubes using a coupled numerical simulation method by comparing the simulation results with analytical solutions. A three-dimensional fluid-structure interaction numerical model is built using Ansys software. Wave propagation is investigated by applying a pressure pulse at the inlet of a fluid-filled elastic tube. The speed of the pressure wave and the radial displacement of the tube are simulated and compared with theoretical values. Simulation results yield a high level of accuracy. Different structural elements are used to represent the tube, and their impact on the results is discussed. The effects of tube material, tube constraints and fluid properties are also investigated in this study.

Introduction

The pulse wave velocity (PWV), the pressure wave propagation speed, in a fluid-filled elastic tube is an important variable in the motion of such tubes. Experimentally the PWV can be determined by using pressure sensors placed at a known distance along a vessel. In blood flow, PWV is considered clinically as a measure of the performance of the arteries’ stiffness, which explains the large amount of work performed in this area. Besides being an excellent test case for FSI simulation, this problem is of direct interest. Our applications include modelling of an apparatus being developed to study digestion ¹ and ducts containing replacement heart valves. ²

The equations relating to PWV have been developed for over a century. The most famous formula is the Moens-Korteweg equation,^3,4 developed in 1878. They found that PWV depends on not only the elastic modulus of the structure, but also the geometrical dimensions of the tube and the fluid density. This equation is derived based on many assumptions and only accounts for the most important factors. A critical review of the equation is given by Lambossy. ⁵ Studies were carried out later by other researchers to provide more precise information on wave propagation characteristics.^6–11 Most studies focused on estimating the elastic properties of the arterial wall and calculating the velocity of the blood flow through the measured PWV. Bergel^12–14 modified the Moens-Korteweg equation by including two more parameters in the calculation, which are the Poisson's ratio and the ratio of thickness to outer radius of the tube. This modification compensates for the errors introduced by the use of the thin wall assumption of the Moens-Korteweg equation.

In the calculation of PWV, the arteries are assumed to be constrained axially, meaning that tube elongation is not considered. Therefore, the effect of tube constraints is not considered when estimating PWV in the circulatory system. However, it has been shown that the wave propagation velocity is also related to the longitudinal support of the tube.^15–19 Different tube constraints can cause different fluid behaviours in an elastic tube. This finding is very important, particularly for the water hammer phenomenon in pipeline systems.

In this study, an investigation of wave propagation in a fluid-filled elastic tube is carried out through numerical simulation. The aim of this study is to verify the accuracy of an FSI model by comparing the simulation results with theoretical values. Firstly, different structural elements are used to build the elastic structure, and their behaviour for various wall thicknesses is discussed. Secondly, the effect of the material properties, which include the wall stiffness and Poisson's ratio, are analysed. Thirdly, some other parameters that are not present in the PWV equations are tested in the numerical model. These parameters include the density of the wall material, the fluid viscosity, and the ramp rate of the pressure pulse. The effect of all parameters is studied for two different longitudinal restraint conditions.

Simulation approach

The validation was performed using Ansys software, in which separate fluids and structural solvers were used iteratively, controlled by Ansys System Coupling to control data transfer between models and convergence. The conservation equations are presented next and then the computational mesh, boundary conditions and model setup.

Fluid model equations

The mass and momentum conservation equations used are shown below for incompressible flow with a moving mesh

$$\begin{array}{c}\mathrm{\nabla }\cdot (\mathit{u}-{\mathit{u}}_g)=0\end{array}$$

(1)

$$\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }(\rho \mathit{u})}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\rho (\mathit{u}-{\mathit{u}}_g)\otimes \mathit{u})=-\mathrm{\nabla }P+\mathrm{\nabla }\cdot \mu (\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^T)}\end{array}$$

(2)

where $\mathit{u}$ is the fluid velocity, ${\mathit{u}}_g$ is the velocity of the moving mesh, $\rho$ is the fluid density, t is time, $\mu$ is the dynamic viscosity, and P is the pressure.

Mechanical model equations

The equation of motion for a linear elastic, homogeneous, isotropic medium is given by:

$$\begin{array}{c}{\rho }_s{\displaystyle \frac{{\mathrm{\partial }}^2\mathit{u}}{\mathrm{\partial }{t}^2}=f+\mathrm{\nabla }\sigma }\end{array}$$

(3)

where ${\rho }_s$ is the instantaneous solid density, $\mathit{u}$ is the displacement vector, $\mathit{f}$ is the body force (henceforth ignored), and $\mathit{\sigma }$ is the stress tensor.

Using conservation of mass, the instantaneous density, ${\rho }_s$ , can be related to the initial material density, ${\rho }_0$ , via

$$\begin{array}{c}{\rho }_s={\rho }_0\det (F)\end{array}$$

(4)

where det denotes the determinant and F is the deformation gradient, which is defined by

$$\begin{array}{c}F={\displaystyle \frac{\mathrm{\partial }\mathit{u}}{\mathrm{\partial }\mathit{X}}}\end{array}$$

(5)

where $\mathit{X}$ is the location of the undeformed model.

In addition, it is readily shown that strain tensor, $\mathit{e}$ , can be calculated from the displacement via

$$\begin{array}{c}e={\displaystyle \frac{1}{2}(\mathrm{\nabla }\mathbf{u}+{(\mathrm{\nabla }\mathbf{u})}^T)}\end{array}$$

(6)

For an isotropic media, the stress tensor is given by

$$\begin{array}{c}\sigma =\lambda \mathrm{tr}(\mathit{e})I+2Ge\end{array}$$

(7)

where I is the identity tensor, tr denotes the trace of a tensor and $\lambda$ and G are the Lamé constants. However, it is usual to work with the Young's modulus, E, and the Poisson's ratio, $\nu$ , that are related to the Lamé constants as follows:

$$\begin{array}{c}\lambda ={\displaystyle \frac{E\nu }{(1+\nu )(1-2\nu )}}\end{array}$$

(8)

$$\begin{array}{c}G={\displaystyle \frac{E}{2(1+\nu )}}\end{array}$$

(9)

where G is the shear modulus (usually written as $\mu$ , but $\mu$ is used for the fluid viscosity in this study).

Using equations (3), (8) and (9), the equation of motion can be rewritten as

$$\begin{array}{c}{\rho }_s{\displaystyle \frac{{\mathrm{\partial }}^2\mathit{u}}{\mathrm{\partial }{t}^2}=(\lambda +2G)\mathrm{\nabla }(\mathrm{\nabla }.u)-G\mathrm{\nabla }\times \mathrm{\nabla }\times u}\end{array}$$

(10)

Setup of the numerical model

The test case considered is the propagation of a pressure wave along a fluid-filled elastic tube that is subjected to a linear pressure ramp at the inlet. The geometry considered is shown in Figure 1. Details of the geometry and material properties are given in Table 1.

Figure 1.

Geometry of the numerical model: showing the fluid region (coloured in blue) and solid region (coloured in orange). The tube is 80 mm long.

Table 1.

Geometrical parameters, fluid and mechanical properties used in the baseline model.

Geometrical dimensions
Radius	$R$	0.002 m
Thickness	$h$	0.0001 m
Length	$L$	0.08 m
Ratio of wall thickness to outer radius	$\gamma$	0.048
Ratio of wall thickness to inner radius	$y$	0.05
Wall properties
Density	${\rho }_s$	1000 kg/m3
Young's modulus	$E$	1 × 106 Pa
Poisson's ratio	$\nu$	0.4999
Fluid properties
Dynamic viscosity	$\mu$	0.001 Pa·s
Density	$\rho$	1000 kg/m3
Pressure ramp rate	$k$	35,250 Pa/s

Fluid computational mesh

The fluid model is developed using ANSYS Fluent, version 2021R1. A size control of 0.5 mm is used to generate the mesh in the fluid domain, shown in Figure 2. It has 14,720 nodes and 18,796 cells with three inflation layers. The minimum orthogonal quality of the mesh is 0.72, and the maximum skewness is 0.28. Mesh refinement tests showed this level of mesh refinement to be sufficient.

Figure 2.

The computational mesh on the fluid domain.

Fluid boundary conditions and setup

The entrance of the tube is set as a pressure inlet with a linearly increasing pressure applied. The exit is set to a static pressure of 0 Pa. The no-slip condition is applied at the tube wall. A dynamic mesh model is enabled, which allows the shape of the domain to change with time. The diffusion-based smoothing method is chosen to diffuse the boundary motion throughout the interior mesh. The inlet and outlet of the tube are allowed to deform in the planes normal to the axial direction. The tube wall is set as a system coupling zone, which transfers data to and from the structural model.

The second-order implicit transient, pressure-based solver is used to conduct the analysis. The coupled scheme is used, which means the momentum and pressure correction equations are solved in a fully coupled manner. Gradients are determined using the least-squares cell-based option, the pressure is determined using the second-order method, and the second-order upwind scheme is used for the momentum equation. Converged solutions are achieved when the residual values for continuity, x, y and z velocities are below 10^–5. The residual values are calculated based on the locally scaled root mean square (RMS) errors.

Mechanical mesh

The choice of meshing strategy is an important component of this study. Of particular interest is an investigation of how the choice of mesh element type affects the solution, particularly in terms of the impact of wall thickness. Three types of elements can be used in Ansys Mechanical to construct the tube: shell, solid, and solid shell. Figure 3 shows the recommended selection criterion, which depends on the aspect ratio, γ, together with a description of the element types.

Figure 3.

Selection criteria and element types available to solve this FSI problem. (a) Element type selection criteria. Based on Akin. ²⁰ (b) Structural element comparison based on the Ansys manual. ²¹

Shell elements are two-dimensional (2D) and are only suitable for modelling thin structures. In a shell model, the structure is represented by a single plate at the mid-surface. Without the thickness dimension, shell elements cannot achieve precise stress and strain profiles through the wall thickness. Each node of the shell element has six degrees of freedom to compensate for the geometric limitation. In addition to the degrees of freedom in translation, the shell element has degrees of freedom in rotations. Therefore, shell elements can capture the bending behaviour at each node explicitly.

Solid elements are well suited for modelling three-dimensional (3D) thick structures because they can simulate the stress and strain profile through the thickness. However, solid elements cannot capture bending behaviour due to the lack of degrees of freedom in rotations. Shear locking occurs when solid elements are used for thin walls. It means that the structure becomes artificially stiffer than its actual stiffness.

Solid shell elements are designed for simulating thin structures. The solid shell model represents the solid structure using a single-layered solid element. The single-layered solid element allows the model to obtain the stress and strain profile through the thickness. Solid shell elements can overcome shear locking when their thickness is small enough due to their unique kinematic formulations. ²²

Different element types are used to build the same model to test their effect. An element size of 0.5 mm is used in all three models. The shell model is represented by a single plate. The solid model has three layers, while the solid shell model only has one. The cross-sectional mesh constructed with different element types are shown in Figure 4 and Table 2 gives details of the number of nodes and elements.

Figure 4.

Cross-sectional mesh on the structural domain for a tube thickness of 0.1 mm (a) shell; (b) solid; (c) solid shell.

Table 2.

Mesh statistics for the different structural models.

	Shell	Solid	Solid shell
Number of nodes in the model	4502	18,032	9016
Number of elements in the model	4480	13,440	4480

Mechanical boundary conditions and setup

The inner surface of the wall is set as a fluid–solid interface, which transfers information between the solid and fluid. Axial displacement is constrained in the structure. Large deflections are enabled to allow accurate calculation of the material deformation. The sparse matrix direct solver is used.

Coupling strategy and setup

Ansys System Coupling, version 2021R1, was used to couple the fluid and structural analyses. A two-way data transfer is performed at the interface. The pressure impulse at the inlet causes the fluid to flow along the tube. The force at the fluid interface is then transferred across the structural interface, causing the structural model to deform. The structural displacement information is then transferred back to the fluid domain. This process is continued until convergence is obtained.

The analysis setup for the FSI model is controlled in the System Coupling setup. The duration of the analysis was set as 0.006 s. A timestep of 0.0005 s was selected after assessing the timestep effect upon results. The maximum number of coupling iterations was set to 20, with typically 13 needed for convergence. Quasi-Newton stabilisation was activated to achieve convergence. This method is derived from an interface quasi-Newton coupling algorithm with an approximation for the inverse of the Jacobian (IQN-ILS). ²³ The initial relaxation factor applied during the start-up iteration was set to 0.01. The maximum number of time steps to be retained in the Quasi-Newton history was set to one. Data transfers at the interface were deemed to have converged when the globally scaled RMS of the residual values are below 0.01.

Ten cores of Intel(R) Xeon Bronze 3204, @1.90 GHz processor with 64 GB RAM were used to run the model. Two cores were used for the Mechanical solver and eight cores were used for the Fluent solver during the simulation.

Analytical solutions

Wave speed

The transient propagation of fluid flow and stress in a thin-walled elastic tube has been studied by many researchers. The most well-known equation is that of Moens-Korteweg^3,4

$$\begin{array}{c}{c}_0=\sqrt{{\displaystyle \frac{Eh}{2\rho R_i}}}\end{array}$$

(11)

where $c_0$ is the wave speed, $\rho$ is the fluid density, $R_i$ is the inner radius of the tube, h is the tube wall thickness, and E is the Young's modulus of the wall material.

The assumptions made in the derivation of the equation are: (i) the tube is filled with incompressible and inviscid fluid, (ii) the wall thickness is very small compared with the radius, (iii) the wall material is elastic and isotropic, (iv) the effect of the Poisson's ratio of the wall material is not considered, and (v) the wall only deforms radially.

Bergel ¹³ modified the Moens-Korteweg equation by including the effect of Poisson's ratio through the wall thickness. This correction compensates for the errors resulting from the thin-wall assumption in the original equation. Although this correction accounts for the effect of the Poisson's ratio, it is still assumed that the length of the tube does not alter, that is, no axial displacement. The modified wave speed becomes

$$\begin{array}{c}{c}_1=c_0\times \sqrt{{\displaystyle \frac{(2-\gamma )}{(2-2\gamma (1-\nu -2{\nu }^2)+{\gamma }^2(1-\nu -2{\nu }^2)-2{\nu }^2)}}}\end{array}$$

(12)

where $\nu$ is the Poisson's ratio, and $\gamma =h/R_o$ , which is the ratio of the wall thickness to the outer radius.

In the derivation of both the Moens-Korteweg equation and Bergel's correction, the elastic tube is not allowed to elongate. The effect of different longitudinal tube constraints was investigated by Wylie et al. ¹⁵ The wave speed formulas are modified based on whether longitudinal motion is allowed or not. Different coefficients are used based on the ratio of the wall thickness to the inner radius. When the ratio $y=h/R_i$ is smaller than 0.08, the tube is considered to have a thin wall. Otherwise, the wall is treated as being thick.

Zero axial displacement. When the tube is anchored against axial movement throughout, the tube is not allowed to elongate. The axial unit strain, ${\xi }_{\mathrm{axial}}=0$ . This scenario is the same as that used to derive equation (12). However, a different approximation is made to that of Bergel, and the following equations are derived.

For a thin-walled tube, the wave speed is

$$\begin{array}{c}{c}_2={\displaystyle \frac{c_0}{\sqrt{1-{\nu }^2}}}\end{array}$$

(13)

For a thick-walled tube, the wave speed is

$$\begin{array}{c}{c}_3={\displaystyle \frac{c_0}{\sqrt{y(1+\nu )+{\displaystyle \frac{2(1-{\nu }^2)}{2+y}}}}}\end{array}$$

(14)

Zero axial stress. When axial movement is allowed along the tube, the axial stress, ${\sigma }_{\mathrm{axial}}=0$ Pa.

For a thin-walled tube, the wave speed is

$$\begin{array}{c}{c}_4=c_0\end{array}$$

(15)

For a thick-walled tube, the wave speed is

$$\begin{array}{c}{c}_5={\displaystyle \frac{c_0}{\sqrt{y(1+\nu )+\left({\displaystyle \frac{2}{2+y}}\right)}}}\end{array}$$

(16)

Radial displacement

In an elastic tube, when a pressure pulse is applied, the radial displacement along the tube can be calculated using the equation below derived from the analytical solution of Womersley.^8,10 However, this equation only applies when the tube is constrained axially.

$$\begin{array}{c}{\delta }_r={\displaystyle \frac{p(z,t)R^2(1-{\nu }^2)}{Eh}}\end{array}$$

(17)

where $p(z,t)$ is the pressure distribution at axial distance z, at time t.

In the simulation, a linearly increasing pressure is applied at the inlet of the fluid domain. The pressure at the inlet follows the following equation:

$$\begin{array}{c}p(0,t)=kt\end{array}$$

(18)

where k is a constant. Therefore, when the tube is constrained axially the radial displacement at the inlet can be calculated from

$$\begin{array}{c}{\delta }_{r0}(t)={\displaystyle \frac{kt{R}^2(1-{\nu }^2)}{Eh}}\end{array}$$

(19)

Effect of fluid compressibility

The fluid is assumed to be incompressible in the Moens-Korteweg equation. By taking the fluid compressibility into account, the wave speed can be calculated using: ⁵

$$\begin{array}{c}{c}_6=\sqrt{{\displaystyle \frac{1}{\rho \left({\displaystyle \frac{1}{K}+{\displaystyle \frac{2R}{Eh}}}\right)}}}\end{array}$$

(20)

The relationship between $c_6$ and $c_0$ is:

$$\begin{array}{c}{c}_6=\sqrt{{\displaystyle \frac{1}{\left({\displaystyle \frac{\rho }{K}+{\displaystyle \frac{1}{c_0}}}\right)}}}\end{array}$$

(21)

where K is the bulk modulus of fluid.

The density of common fluids is of the order of 10³ kg/m³, and the bulk modulus is of the order of 10⁹ kPa. Therefore, $\rho /K$ is very small compared with the inverse propagation speed, and it is reasonable to assume the fluid to be incompressible when calculating the wave speed. However, it should be noted that if gas is present in the fluid, the effect of the mixture compressibility becomes significant, and the wave speed can be reduced dramatically. Detailed calculations can be found in Kobori et al. ²⁴ and Pearsall. ¹⁶

Model verification

Effect of structural element type and wall thickness

The effect of the structural element type is examined by building the same model using different element types. Figure 5 shows the radial mesh displacement along the tube for shell elements. Figure 5(a) shows that the pressure impulse at the inlet increases over time, causing the wall displacement to gradually increase. The simulation is terminated before the wave reaches the end of the tube as it would cause reflections and disturb the profiles. The radial displacement at the inlet over time is recorded and compared with the theoretical results following equation (19). The results show a high level of accuracy for different wall thicknesses.

Figure 5.

Propagation behaviour using shell elements: (a) Radial mesh displacements along the tube, (b) Radial displacements at the inlet compared with the analytical solution for various wall thicknesses.

The wave speeds are estimated by calculating the distance travelled by the wave over time. The calculated wave speeds for different structural models are shown in Figure 6. The simulation results from all models are very similar when $y$ is smaller than 0.1. The trends begin to deviate as the wall gets thicker. When $y>0.1$ , $c_3$ provides a better alignment with the solid and solid shell models compared with the shell model. These results are consistent with the structural element selection criteria shown in Figure 3. Solid shell and solid elements provide more accurate results when the wall is relatively thick. Moreover, the solid shell model and solid model produce very similar results. Therefore, the solid shell model is used to represent a solid model in the subsequent studies because it is much less computationally expensive.

Figure 6.

The effect of wall thickness on wave speed for different element types.

Mesh and timestep independence studies

A mesh-resolution study was performed by reducing the element size from 0.5 to 0.25 mm. Simulation results are compared against the analytical propagation speed $c_1$ as it provides a better result when y is smaller than 0.1. The results showed that both models give the same solution, demonstrating mesh independence in both models. ¹

A timestep independence study in which the timestep was reduced to 0.2 from 0.5 ms also showed the original result to be timestep independent. ¹

Effect of wall stiffness

The effect of wall stiffness is examined in the simulations and the results for the propagation speed are compared with the analytical value $c_1$ for a case in which the tube is constrained axially. As shown in Figure 7, both shells and solid shells show good agreement. The wave speed increases as the wall stiffness increases and the data fit perfectly with equation (12). The radial displacement at the inlet for different wall stiffnesses fits very closely with equation (19).

Figure 7.

The effect of wall stiffness when the tube is constrained axially on: (a) wave speed; (b) radial displacement.

When the tube is constrained axially, axial strain along the tube is zero, while the axial stress varies along the tube. The stress profile matches very well with the assumption of Case (a) from Wylie et al. ¹⁵ (see ¹ ).

When the tube is free to elongate, the simulated wave speeds are compared with the values for wave speed $c_4$ and are shown in Figure 8. Without the axial constraint, the wave speeds become slower, while the radial displacement becomes larger. The axial stress along the tube is zero, while the axial strain varies along the tube. These profiles also match well with the assumption of Case (b) from Wylie et al. ¹⁵ (see ¹ ).

Figure 8.

The effect of wall stiffness on wave speed when the tube is free to elongate on: (a) wave speed; (b) radial displacement.

Effect of Poisson's ratio

The effect of Poisson's ratio is examined when the tube is constrained axially. Figure 9 shows that as the Poisson's ratio increases, wave speeds from both element types increase, following the trend of the analytic solution yielding the wave speed $c_1$ . When the Poisson's ratio increases from 0 to 0.5, the wave speed increases by 15% in the shell model and 6% in the solid shell model. The radial displacement at the inlet for different wall stiffness matches well with equation (19).

Figure 9.

The effect of Poisson's ratio when the tube is constrained axially on: (a) wave speed; (b) radial displacement.

When the tube is free to elongate, the Poisson's ratio does not have a significant impact on either the wave speed or the radial displacement. When the Poisson's ratio is increased from 0 to 0.5, the wave speed for both models decreased by only 2%, as can be seen in Figure 10.

Figure 10.

The effect of Poisson's ratio when the tube is free to elongate on: (a) wave speed; (b) radial displacement.

Effect of wall material density

In the theoretical PWV calculation, the effect of the added mass, which is the wall mass and the fluid carried with it, is not considered. In the study of Morgan and Ferrante, ⁷ the effect of solid density was discussed. The study was then extended by Womersley, ¹⁰ who examined the effect of mass-loading of the solid. He stated that the wall material density only affects the wave speed if the tube constraint allows longitudinal motion. Womersley ¹⁰ defined the term $K_{ml}$ to include the effect of mass-loading and longitudinal constraint:

$$\begin{array}{c}{K}_{ml}=\left(1+{\displaystyle \frac{h_s}{h}\cdot {\displaystyle \frac{{\rho }_s{R}_s}{\rho R_i}}}\right)\left(1-{\displaystyle \frac{m^2}{n^2}}\right)\end{array}$$

(22)

where $K_{ml}$ is related to the wave speed, $h_s$ , ${\rho }_s$ and $R_s$ are the thickness, density and mean radius of a tube that provides the added mass of the solid, m is the natural frequency of the longitudinal constraint and n is the frequency of the pressure pulse.

If $m=n$ , $K_{ml}=0$ , the tube behaves as if there is no mass loading. If the longitudinal constraint is stiff enough to prevent any longitudinal motion of the wall, $m\gg n$ and $K_{ml}\to -\mathrm{\infty }$ . In this case, the value of ${\rho }_s$ does not have any effect. However, if the constraint is not very stiff and leads to longitudinal motion, $K_{ml}$ is related to the ratio of ${\rho }_s/\rho$ . In this section, the effect of the wall material density is estimated through numerical simulations.

When the tube is constrained axially, results for the effect of the wall material density are shown in Figure 11. The largest difference between the results is only 3%. Therefore, the wall material density does not have a significant impact on either the wave speed or radial displacement when the tube is constrained axially.

Figure 11.

The effect of wall material density when the tube is constrained axially on: (a) wave speed; (b) radial displacement.

When the tube is free to elongate, the wall material density shows an effect on both wave speed and radial displacement, as shown in Figure 12. When the density increases from 1000 to 10,000 kg/m³, the wave speed decreases by around 16%. It also has an impact on the radial displacement. However, the radial displacement at the inlet is almost the same when the density increases from 5000 to 10,000 kg/m³.

Figure 12.

The effect of wall material density when the tube is free to elongate on: (a) wave speed; (b) radial displacement.

The effect of fluid viscosity

The fluid is assumed to be inviscid in the PWV calculation to simplify the equations. However, the fluid viscosity affects the wave velocity significantly according to Womersley ⁸ and Bergel. ¹² In this section, the effect of the fluid viscosity is examined through numerical simulation. As shown in Figure 13, the fluid viscosity has a significant impact on wave speed when the tube is constrained axially. There is a 10% difference when the viscosity is changed from 0.001 to 0.1 Pa·s. However, the radial displacement is independent of the fluid viscosity.

Figure 13.

The effect of fluid viscosity when the tube is constrained axially on: (a) wave speed; (b) radial displacement.

The wave speed is also affected by fluid viscosity when the tube is free to elongate. When the fluid viscosity is increased from 0.001 to 0.1 Pa·s, the wave speed decreases by around 16% as can be seen in Figure 14. Also, the radial displacement decreases when the fluid viscosity increases up to 0.01 Pa.s, and then the effect is small.

Figure 14.

The effect of fluid viscosity when the tube is free to elongate on: (a) wave speed; (b) radial displacement.

The effect of pressure ramp rate

The effect of the pressure ramp rate is examined when the tube is constrained axially. As shown in Figure 15, the wave speeds remain the same in both models when the pressure ramp rate changes. The radial displacement results fit perfectly with results from equation (19). Similar behaviour is observed for an unconstrained tube.

Figure 15.

The effect of pressure ramp rate when the tube is constrained axially on: (a) wave speed; (b) radial displacement.

Conclusions

Model verification is crucial for simulation development, yet it is often neglected, especially when commercial software is used. In this study, the accuracy of the Ansys System Coupling approach to model FSI of wave propagation in a fluid-filled tube is verified. The choice of structural element type is shown to be important when modelling different wall thicknesses and for moderately thick walls solid shells gave the best results. Specifically, if the ratio of the wall thickness to the tube inner radius is less than 0.1 (y < 0.1), shell elements perform best and above this ratio solid shell or solid elements are required. This corresponds closely with the limit set by Wylie et al. ¹⁵ of y < 0.08 for a tube to have thin walls. Where applicable, solid shells are preferred over solids due to their lower computational cost.

As would be expected the tube constraint has an impact on axial stress and strain profiles of the elastic tube material, and this manifests itself in changed wave speeds and wall displacement for different tube materials. This is an important observation, as in many real systems the choice of boundary condition is not obvious. In biomedical systems, there are very few cases where a simple fixed displacement condition can be applied. It is important to examine the sensitivity to both the location of the constraint and the type of boundary condition. This ensures the choice leads to physical results and determines the sensitivity to the condition applied.

Key parameters such as wall thickness and material properties were varied, and results showed excellent agreement with analytical predictions. The simulation results also show the effect of the tube materials, tube constraints and fluid properties, which are not readily available from theoretical studies.

Importantly, it can be concluded that the solution methodology presented here is suitable for the simulation of peristaltic flows that are important in a wide variety of biomedical and engineering simulations. The validation presented here is for a straight tube and demonstrates the accuracy of the FSI coupling methodology used. The same algorithm and methodology can be applied to scans of PWV or other flexible tubes and the same accuracy for properly resolved simulations can be expected. Previous work ²⁵ has shown that the fluid model is very accurate for simulations when a known travelling wave motion is applied to the wall, so the combination of this work and the previous work validates the modelling approach for a wide range of conditions.