In this manuscript, we model the transmission of pressure pulses along a straight pipe made of a viscoelastic material, filled with a fluid. This model can be separated in two parts: The first one refers to the dynamic stress/strain behavior of the material of the pipe, also known as material model, and the other one refers to the wave equation governing the pressure wave transmission. We have compared this model to a set of experimental strain/stress and strain/time curves and adjusted the parameters for a polymer for medical devices, polydimethylsiloxane.
We have used two well-known viscoelastic material model, the Maxwell model (Fig. 1a) and the Zener model (Fig. 1b). The Maxwell model is defined by Eq. (1), where is the strain and the stress. and are the static Young modulus and the viscous factor, respectively. The Zener model (also known as standard linear viscoelastic model or three-parameter model) is defined by Eq. (2), where .
The hypotheses assumed for the model are that the fluid is incompressible; the effect of the viscosity of the fluid is negligible; and the deformation of the tube is small related to its radius. Fluid movement is a plug flow, so the velocity does not depend on the radial dimension. In fact, a published study in which speed was locally measured (Pezzinga et al. 2014) shows fairly flat velocity profiles in the hydraulic transients. However, it should be noted that this study was carried out for turbulent regimes, whereas flow is usually laminar in the circulatory system.
For a pipe with wall thickness (e) small compared to the inner radius (r), the relationship between the pressure difference between the pressure inside and outside the pipe (p) and the tangential stress () at the pipe wall (Fig. 1c) is:
The continuity equation for a differential volume inside the tube gives:
In addition, the balance of the momentum gives:
where u is the velocity of the fluid and its density. In the last term, fis the frictional force per unit length that depends on the viscosity and the shape of the velocity (Willemet and Alastruey 2015). The term of convective acceleration may be negligible when the velocity is low. In this work, the convective acceleration and frictional force terms have been neglected.
If the strain is , combining Eq. (1) with Eq. (4) gives:
Therefore, Eqs. (5, 6), after cross-derivation, give the unique equation:
Taking into account (3), if c is the pressure wave velocity, the model results in:
This is a wave equation with a linear dissipation term. When tends to infinity, does too and Eq. (8b) becomes the well-known equation for (undamped) pressure pulses inside elastic pipes, like in a water hammer.
When the initial and boundary conditions are convenient, Eq. (8b) can be solved by separation of variables and Sturm–Liouville series. For the conditions set in (9), this becomes an eigenvalue problem that allows calculating the decaying rate of the pressure wave. These eigenvalues define the inherent damping ability of the material.
The eigenvalues, which are the spatial frequencies, are: , and the formal solution to (8) is:
where with the discriminant
Note that are always negative. Hence, if the discriminant is negative, the corresponding term is underdamped. If the discriminant is zero, the corresponding term is critically damped. If the discriminant is positive, the corresponding term is overdamped. Note that the sequence of the values of this discriminant (11) is decreasing and unbounded. In consequence, there will be infinite negative terms (underdamped). However, for low values of n, the discriminant could be positive and hence the lowest spatial frequencies may be overdamped. This is when
For underdamped cases, the solution to (8) can be written as a function of time and space:
where is the position of the pressure pulse and is the spatial coordinate.
If and are constants to be determined from the initial conditions, the auxiliary function is defined as follows:
is an undamped wave, and most importantly a periodic function in both time and space. Consequently, it does not show any irreversible decay, and therefore, the decay rate of the wave is the negative exponential in (13a) with constant of time .
For those cases where condition (12) is not fulfilled, a number of overdamped terms arise. This happens for values of n lower than a critical value k, which is . Then, the solution becomes:
Note that , so the wave will show an overdamped part and an underdamped part. One part of the overdamped terms decay faster than the undamped terms, and the other part decays slower. The dominant decay rate of the overdamped part is controlled by .
Time derivation of the Zener model (2) gives:
The continuity Eq. (4) gives:
And, after time derivation:
Combining (15) with (16, 17), one has:
And, deriving respect to time once again:
The balance of the momentum, assuming no convective acceleration and no friction (5) and deriving respect to space gives:
And deriving again respect to time:
If and , one can now combine Eq. (19) with (3, 20 and 21):
This third-order partial differential equation appears as a linearized model in nonlinear acoustics, under the name of Moore–Gibson–Thompson equation. Taking the dominant terms (those of third order), this is a wave equation with a wave velocity c. We propose the next formal solution, which is analogous to the one found in (13):
where is the solution of the equation:
and
When pressure is zero, (25) gives and the function is:
where L is the length of the spatial domain and every constant. From (25) one gets the cubic characteristic equation, for the unknown :
This polynomial does not have any negative real root. Using the formulae of Cardano, the roots of this polynomial are:
where
The lowest value of the real part of controls the decaying rate and that is:
The value of decreases as n increases and tends to a constant value depending on the parameters of the material model:
Therefore, the limiting value of the constant can be stated from (30) and (31):
This limiting value is known as the “essential spectrum” of the problem. Observe in the expression (30) that, for a range of cases, the value of is almost independent on n and the exponential part of the solution (23) can be taken out of the sum. This limiting value of can be approximated by (32), and this is when:
The same is true when the pulse width is very short compared with the length of the space considered, because the lowest harmonic terms in (23) are very weak.
Note that the imaginary part of (28b,28c) represents the harmonics of an oscillation with increasing frequency. Therefore, the solution for (23) can be written as a function of time and space:
where is the position of the pressure pulse and is the spatial coordinate. is a bounded function which is periodic in both time and space or damped in time and periodic in space, and can be calculated by (35).
When condition (33) is not fulfilled, the value of in Eq. (30) depends on the value of n. In that case, the exponential part of the solution (23) can’t be taken out of the sum and the decaying rate becomes a weighted sum of n exponential decays limited by (lowest decaying rate) and (highest decaying rate).
In summary, expressions (13) and (34-35) give the decaying rate of a pressure wave along a tube when condition (33) is fulfilled only for the Zener model. Note that both expressions refer to the pressure, which is proportional to stress (3); therefore, the decaying rate for the stress is the same for pressure and stress. When (33) is not fulfilled, the decaying rate is bounded between the decaying rate calculated for and the one calculated for with Eq. (30).
The formulae previously developed are compared to numerical one-dimensional simulations. We propose the problem with a realistic blood pressure pulse arising from a heartbeat. The pressure wave propagates along a straight viscoelastic pipe of length L (Fig. 1d), with pulsatile inlet f(t) as shown in Eq. (36):
In consequence, for the Maxwell material model, the set of conditions of the problem are defined in (37) as:
The partial differential equation (PDE) is second order both in time and in space. It has been solved by the method of lines (Saucez et al. 2004). The space dimension has been discretized using finite differences, by centered formulas of second order. This gives a system of ordinary differential equation (ODE) in time, where h is the space step of the discretization.
To reduce these second-order ODEs to first-order ones, one takes:
This gives 2N equations:
And the conditions:
The ODE system (40,41) is solved using the Dormand–Prince method (Dormand and Prince 1980), programmed in MATLAB (ode45).
For the Zener model, the set of conditions of the problem are defined in (42), analogous to (37), with the pressure pulse defined previously in (36):
Taking into account (39) and (42):
The PDE gives the system of 3N ODE ():
And the boundary conditions give:
The ODE system (44,45) is solved using the same MATLAB code.
The ability to damp pressure waves was tested in polydimethylsiloxane (PDMS, Dow Corning, Midland, MI), a common material used in medical applications. Maxwell and Zener models were adjusted to represent the viscoelastic behavior of PDMS using a stress/strain experiment. Stress/strain tests were carried out using Instron ElectroPulsTM E3000a tensile rig. Force was obtained from load cells and displacement obtained from grip displacement. Briefly, a uniaxial sinusoidal displacement at 1Hz was applied to a flattened half cylinder of PDMS of 14.28 mm diameter and 2.37 mm thickness. This tube was cut into slices of 10 mm width and straightened to a length of 37.38 mm. The PDMS sample was held to the uniaxial test machine allowing a length between fixations of 20 mm. (See supplemental Video 1) Tube flattening provoked a 20% compression on the exterior face and a 20% tension on the interior face. The displacement applied deformed the sample up to 7% to obtain stress/time and strain/time curves. Our experiment measured the average tensile strength along the tube. We performed a series of simulations to confirm that the significant thickness of our pipe (0.166 thickness-to-diameter ratio) would not affect our results. Indeed, we observed that, for ratios between 0.1 and 0.2, the error between theoretical and average stress ranged approximately .
The stress/strain experiment was simulated in parallel to estimate the material models’ parameters that fit best the experimental data. These simulations were run using the Dormand–Prince method, with the experimental strain/time curve as an input and the stress/time curves as an output. The error between the experimental and the simulated stress/time curves was calculated as the average of the orthogonal distance from each point of the experimental data to the simulated curve. This error was minimized using the Nelder–Mead simplex method (Lagarias et al. 1998) to find the material models’ parameters that best represented the material behavior.
