Mathematical problem formulation
Derivation of 1D hyperbolic system
The fluid mechanics model that is developed and solved in this package is based on the equations of mass and momentum with density, temperature and viscosity considered as constant. Moreover, it is assumed that all the elements of the model lie on the same height; as a result, the gravitational forces are cancelled. Lastly, the flow is considered to be axisymmetric.
For thin cylinders and axisymmetric flow, the Navier Stokes equations are expressed as follow:
- r direction
- x direction
where \( v \) and \( u \) is the velocity of the fluid along r (radius) and x (cylinder axis), respectively. Thus, \( v \) and \( u \) depend on \( r \), \( x \) and \( t \) (time) alternatively expressed as \( v = v(r, x, t) \) and \( u = u(r, x, t) \). Lastly, \( p \) is the pressure of the fluid, explicitly expressed as \( p = p(r, x, t) \). \( \rho \) and \( \nu \) is the density and the dynamic viscosity of the fluid, respectively.
If Eq. (2) is multiplied with the density, \( \rho \), the following formula is obtained $$ \begin{align} \tag{3} \rho \frac{\partial u}{\partial t} + \rho \left( v \frac{\partial u}{\partial r} + u \frac{\partial u}{\partial x} \right) = -\frac{\partial p}{\partial x} + \mu \frac{1}{r} \frac{\partial}{\partial r} \left ( r \frac{\partial u}{\partial r}\right) \end{align} $$
The equation of mass, for incompressible fluid and axisymmetric flow is expressed as $$ \begin{align} \tag{4} \frac{1}{r} \frac{\partial (rv)}{\partial r} + \frac{\partial u}{\partial x} = 0 \end{align} $$
Eqs. (3) and (4), depend on r. In order to simplify the afforementioned equations, integration along the cylindrical tube radius and multiplication by \( 2 \pi r dr \) is carried out. Therefore, Eq. (3) is transformed as $$ \begin{align} \underbrace{\rho \int_0^R \frac{\partial u}{\partial t} 2 \pi rdr}_{I} + \underbrace{\rho\int_0^R v\frac{\partial u}{\partial t} 2 \pi r dr}_{II} + \underbrace{\rho\int_0^R u \frac{\partial u}{\partial x} 2 \pi rdr}_{III} = \underbrace{-\int_0^R \frac{\partial p}{\partial x} 2\pi rdr}_{IV} + \underbrace{2\pi \mu \left[r \frac{\partial u}{\partial r} \right]_0^R}_{V} \tag{5} \end{align} $$
It well known that in thin-long cylinders with low curvature and axisymmetric flow, the wall shear stresses can be given by the following relationship $$ \begin{align} \tag{6} \tau_R = \mu \frac{\partial u}{\partial r} \bigg\rvert_R \end{align} $$
Substituting Eq. (6) in the term \( V \) of Eq. (5), the following formula is obtained $$ \begin{align} \underbrace{\rho \int_0^R \frac{\partial u}{\partial t} 2 \pi rdr}_{I} + \underbrace{\rho\int_0^R v\frac{\partial u}{\partial t} 2 \pi r dr}_{II} + \underbrace{\rho\int_0^R u \frac{\partial u}{\partial x} 2 \pi rdr}_{III} = \underbrace{-\int_0^R \frac{\partial p}{\partial x} 2\pi rdr}_{IV} + \underbrace{2\pi R \tau_R}_{V} \tag{7} \end{align} $$
To simplify Eq. (7), the following steps can be carried out:
- Term \( \bf{I} \) By using Leibniz integral rule, term \( \bf{I} \) transforms to
The average velocity in a cross-section of the cylindrical tube can be expressed as $$ \begin{align} \tag{9} \bar{u} = \frac{1}{A}\int_0^R u 2 \pi rdr = \frac{Q}{A} \end{align} $$
So, the first term of Eq. (7) can be expressed as (also with \( A=\pi R^2 \)) $$ \begin{align} \tag{10} \underbrace{\rho \int_0^R \frac{\partial u}{\partial t} 2 \pi rdr}_{I} = \rho \frac{\partial (A \bar{u})}{\partial t} - \rho u(R, x, t)\frac{\partial A}{\partial t} \end{align} $$
Where the product \( \bar{u} A \) is the volumetric flow rate \( Q(x, t) \) in a cross-section.
- Terms \( \bf{II} \) and \( \bf{III} \)
and by Eq. (4), the above becomes $$ \begin{align} \rho \int_0^R v \frac{\partial u}{\partial r} 2 \pi rdr + \rho \int_0^R u \frac{\partial u}{\partial x} 2 \pi r dr = 2 \pi \rho \left[ rvu \right]_0^R + 2 \pi \rho \left( \int_0^R ru \frac{\partial u}{\partial x} dr + \int_0^R ru \frac{\partial u}{\partial x} dr \right) \tag{12} \end{align} $$
or alternatively $$ \begin{align} \tag{13} \rho \int_0^R v \frac{\partial u}{\partial r} 2 \pi rdr + \rho \int_0^R u \frac{\partial u}{\partial x} 2 \pi r dr = 2 \pi \rho \left[ rvu \right]_0^R + 2 \pi \rho \underbrace{ \int_0^R \frac{\partial (ru^2)}{\partial x} dr}_{i^{\star}} \end{align} $$
where the last integral \( i^{\star} \) of the last relationship with Leibnitz integral rule becomes $$ \begin{align} \tag{14} \int_0^R \frac{\partial (ru^2)}{\partial x} dr = \frac{\partial}{\partial x} \int_0^R ru^2 dr - \frac{1}{2}u^2(R,x,t) \frac{\partial R^2}{\partial x} \end{align} $$
Substituting Eq (14) to Eq. (13) (and by \( A = \pi R^2 \)) the sum takes the final form $$ \begin{align} \underbrace{\rho \int_0^R v \frac{\partial u}{\partial r} 2 \pi rdr}_{II} + \underbrace{\rho \int_0^R u \frac{\partial u}{\partial x} 2 \pi r dr}_{III} = 2 \pi \rho R v(R, x, t) u(R, x, t) + 2 \pi \rho \frac{\partial }{\partial x} \int_0^R ru^2 dr - \rho u^2(R, x, t)\frac{\partial A}{\partial x} \tag{15} \end{align} $$
Furthermore, integrating Eq. (4) from \( 0 \) to \( R(x, t) \) and multiplying by \( 2 \pi r dr \), continuity align transforms to $$ \begin{align} \tag{16} 2 \pi R v(R, x, t) + \underbrace{\int_0^R \frac{\partial u}{\partial x} 2 \pi r dr}_{c^{\star}} = 0 \end{align} $$
For the integral \( c^{\star} \) of Eq. (16), by using Leibnitz integral rule, the above relationship becomes $$ \begin{align} \tag{17} 2 \pi R v(R, x, t) + \frac{\partial}{\partial x} \int_0^R u 2 \pi rdr - u(R, x, t)\frac{\partial (\pi R^2)}{\partial x} = 0 \end{align} $$
Substituting \( A = \pi R^2 \) and \( \bar{u} = \frac{1}{A}\int_0^R u 2 \pi rdr = \frac{Q}{A} \), Eq. (17) becomes $$ \begin{align} \tag{18} 2 \pi R v(R, x, t) = - \frac{\partial (A\bar{u})}{\partial x} + u(R, x, t)\frac{\partial A}{\partial x} \end{align} $$
Finally, substituting Eq. (18) to the first member of the right part of Eq. (15), the final sum transforms to $$ \begin{align} \underbrace{\rho \int_0^R v \frac{\partial u}{\partial r} 2 \pi rdr}_{II} + \underbrace{\rho \int_0^R u \frac{\partial u}{\partial x} 2 \pi r dr}_{III} = - \rho u(R, x, t) \frac{\partial (A\bar{u})}{\partial x} + 2 \pi \rho \frac{\partial }{\partial x} \int_0^R ru^2 dr \tag{19} \end{align} $$
- Term \( \bf{IV} \)
Taking into account the relationship \( A = \pi R^2 \) along with the average pressure in a cross-section; \( \bar{p} = \frac{1}{A} \int_0^R p 2 \pi r dr = \frac{F}{A} \), the above integral transforms to $$ \begin{align} 2 \pi \int_0^R \frac{\partial (rp)}{\partial x} dr = \frac{\partial (A \bar{p})}{\partial x} - p(R, x, t)\frac{\partial A}{\partial x} \tag{21} \end{align} $$
which by integrating by parts the first term of the right part it transforms to $$ \begin{align} 2 \pi \int_0^R \frac{\partial (rp)}{\partial x} dr = A \frac{\partial \bar{p}}{\partial x} + \bar{p}\frac{\partial A}{\partial x} - p(R, x, t)\frac{\partial A}{\partial x} \tag{22} \end{align} $$
If it is assumed that in the above equation the gradients of pressure, \( p(R, x, t) \), along the R direction are negligible (negligible curvature of the tube), the last two terms are cancelled; as a result, the last equation becomes $$ \begin{align} 2 \pi \int_0^R \frac{\partial (rp)}{\partial x} dr = A \frac{\partial \bar{p}}{\partial x} \tag{23} \end{align} $$
or $$ \begin{align} \tag{24} \underbrace{\int_0^R \frac{\partial p}{\partial x} 2\pi r dr}_{IV} = A \frac{\partial \bar{p}}{\partial x} \end{align} $$
Now, substituting the integral terms of Eq. (7) with the integrals of Eqs. (10), (19) and (24), the x momentum equation transforms to $$ \begin{align} \underbrace{\rho \frac{\partial (A \bar{u})}{\partial t} - \rho u(R, x, t)\frac{\partial A}{\partial t}}_{I} - \underbrace{ \rho u(R, x, t) \frac{\partial (A\bar{u})}{\partial x} + \rho \frac{\partial }{\partial x} \int_0^R u^2 2 \pi r dr}_{II + III} = - \underbrace{A \frac{\partial \bar{p}}{\partial x}}_{IV} + \underbrace{2\pi R \tau_R}_{V} \tag{25} \end{align} $$
If we assume \( \bar{u} = u(R, x, t) \), a new correction factor, \( \beta \), is added to the above momentum equation $$ \begin{align} \tag{26} \beta = \frac{1}{\bar{u}^2 A} \int_0^R u^2 2 \pi r dr \Rightarrow \int_0^R u^2 2 \pi r dr = \beta \bar{u}^2 A \end{align} $$
Thus, substituting Eq. (26) to Eq. (25), the x-momentum equation can be expressed as $$ \begin{align} \tag{27} \rho \frac{\partial (A \bar{u})}{\partial t} - \rho u(R, x, t)\left[ \frac{\partial A}{\partial t} + \frac{\partial (A \bar{u})}{\partial x} \right] + \rho \frac{\partial (\beta \bar{u}^2 A)}{\partial x} = - \frac{\partial \bar{p}}{\partial x}A + 2 \pi R \tau_R \end{align} $$
with \( 1 \leq \beta < \frac{4}{3} \). Normally, \( \beta \) equal to unity and \( \frac{4}{3} \) occurs in turbulent flow (uniform velocity profile) and laminar flow (parabolic profile), respectively. Since the range of \( \beta \) is very narrow along with its presence in convective acceleration term (it is evident that it is has low contribution in contrast to the rest terms of the momentum equation), usually \( \beta \) is regarded to unity. The last approximation only applies when the boundary layer is thin compared to vessel radius. Moreover, this assumption has been proven to be good estimate for pulsatile flows of high Womersley number.
The boundary layer can be described by \( r = R(x, t) \), or alternatively, \( h=r - R(x, t) = 0 \). The total time derivative \( \frac{Dh}{Dt} = 0 \) since the walls are impermeable. Analysing the material derivative we get $$ \begin{align} \frac{\partial h}{\partial t} + v(R, x, t)\frac{\partial h}{\partial r} + u(R, x, t)\frac{\partial h}{\partial x} = 0 \tag{28} \end{align} $$
Since \( \frac{\partial r}{\partial t} = 0 \), the above equation tranforms to $$ \begin{align} \frac{\partial R(x, t)}{\partial t} + v(R, x, t) - u(R, x, t)\frac{\partial R(x, t)}{\partial x} = 0 \tag{29} \end{align} $$
Solving with respect to radial velocity and multiplying both sides by \( 2 \pi R \) (and \( A = \pi R^2 \)) the following equation can be obtained $$ \begin{align} 2 \pi R v(R, x, t) = \frac{\partial A}{\partial t} + u(R, x, t)\frac{\partial A}{\partial x} \tag{30} \end{align} $$
From the last equation and Eq. (18) $$ \begin{align} \frac{\partial A}{\partial t} = - \frac{\partial (A\bar{u})}{\partial x} \tag{31} \end{align} $$
Thus, the 1D continuity equation becomes $$ \begin{align} \tag{32} \frac{\partial A}{\partial t} + \frac{\partial \overbrace{(A \bar{u})}^{Q}}{\partial x} = 0 \end{align} $$
Finally, from the 1D continuity equation, Eq. (32), the x-momentum equation is expressed as $$ \begin{align} \tag{33} \frac{\partial Q}{\partial t} + \frac{\partial (\beta Q^2 / A)}{\partial x} = - \frac{A}{\rho} \frac{\partial \bar{p}}{\partial x} + \frac{2 \pi R \tau_R}{\rho} \end{align} $$
Velocity profile
Assuming that the velocity profile is parabolic, i.e. $$ \begin{align} u_z = 2 u \left( 1 - \frac{r^2}{R^2} \right) \tag{34} \end{align} $$
the last term on the right part of Eq. (33) becomes $$ \begin{align} 2 \pi R \frac{ \overbrace{\mu \left[ \frac{\partial u_z}{\partial r}\right]_R}^{\tau_R} }{\rho} = - \frac{8 \pi \mu q}{\rho A} \tag{35} \end{align} $$
and the correction factor, \( \beta \), in momentum equation becomes $$ \begin{align} \beta = \frac{1}{u^2 A} \int_0^R 4 u^2 \left( 1 - \frac{r^2}{R^2} \right)^2 2 \pi r dr = \frac{4}{3} \tag{36} \end{align} $$
However, many studies ([1]) have shown that the parabolic velocity profile is valid for steady and laminar flow. On the other hand, the system is more dynamic and the velocity profile changes according to the flow conditions. In general, for laminar flow in tappering vessels, the velocity profile is rather flat ([2]). Therefore, a better approach is to assume that the velocity profile is flat with a boundary layer; \( \delta \), described as $$ u_z = \begin{cases} u, & \text{for $r \leq R - \delta$} \\ u(R - r)/ \delta, & \text{for $R-\delta < r \leq R$} \end{cases} $$
In this case, the last term on the right part of Eq. (33) can be written as $$ \begin{align} 2 \pi R \frac{ \overbrace{\mu \left[ \frac{\partial u_z}{\partial r}\right]_R}^{\tau_R} }{\rho} = - \frac{2 \pi \mu u R}{\rho \delta} = - \frac{2 \pi \nu q R}{\delta A} \tag{37} \end{align} $$
and the correction factor, \( \beta \), becomes $$ \begin{align} \beta = \frac{1}{u^2 A} \int_0^{R - \delta} u^2 2 \pi r dr + \frac{1}{u^2 A} \int_{R - \delta}^R \left( \frac{u(R-r)}{\delta} \right)^2 2 \pi r dr \approx 1 \tag{38} \end{align} $$
Therefore, the last approximation applies if the boundary layer is thin compared to vessel radius. This value is often used since it leads to a considerable simplification in the analysis and the loss of relevance of the model is very small in most cases [3].
State equation
As can be observed in the above system, there are three dependent variables, \( A \), \( Q \) and \( p \). Therefore, in order to solve the above system, pressure \( p \) has to be expressed with respect to area \( A \). At this stage, we use a linear elastic model which balances the internal and external forces in radial direction of a surface element on the vessel wall. The pressure-area relation is expressed as $$ \begin{align} \tag{39} p = \frac{4}{3}\frac{\sqrt{\pi}E(x)h}{A_0}\left( \sqrt{A} - \sqrt{A_0} \right) \end{align} $$
with \( E(x) \) and \( h \), the Elastic modulus and the thickness of the vessel wall, respectively.
Olufsen [4] (pp 38) used volume compliance clinical data reported earlier in Stergiopoulos et al. [5] and calculated \( (Eh)/r_0 \) values (from \( (Eh)/r_0 = (3 A_0 L)/(2 C_v) \) with \( C_v \), the volume compliance of a vessel). Thereafter, an exponential empirical model was fitted, expressed as $$ \begin{align} \frac{Eh}{r_0} = k_1\exp{k_2 r_0} + k_3 \tag{40} \end{align} $$
Therefore, we can express the elastic modulus with respect to reference diameter \( r_0 \) via this empirical relationship. The \( \bf{k} \) vector with the structural parameters can then be obtained with a least square optimisation fit. Olufsen run the optimisation and found that the optimised \( \bf{k} \) vector is as follows $$ \mathbf{k} = \left[ \begin{matrix} k_1 \\ k_2 \\ k_3 \end{matrix} \right] = \left[ \begin{matrix} 2.00 \\ -2.253 \\ 0.0865 \end{matrix} \right]\left( \begin{matrix} tonnes/(s^2 mm) \\ mm^1 \\ tonnes/(s^2 mm) \end{matrix}\right) $$
The above empirical relationship can now be inserted in the pressure-are relationship, trasforming Eq (39) as $$ \begin{align} \tag{41} p = f(r_0, \mathbf{k})\left( \sqrt{\frac{A}{A_0}} - 1 \right) \end{align} $$
with \( f(r_0, \mathbf{k}) = \frac{4}{3}\left(k_1\exp{k_2 r_0} + k_3 \right) \).
Conservational form and coupling of PDEs
After the assumption of a flat velocity profile, the system of PDEs that governs the propagation of pulse wave in large arteries is as follows $$ \begin{aligned} \frac{\partial A}{\partial t} + \frac{\partial Q}{\partial x} = 0 \\ \frac{\partial Q}{\partial t} + \frac{\partial (Q^2 / A)}{\partial x} + \frac{A}{\rho} \frac{\partial p}{\partial x} = - \frac{2 \pi \nu q R}{\delta A} \end{aligned} \tag{42} $$
To solve this system numerically, we have to transform it to conservation form. The derivative of pressure with respect to x axis in the momentum equation can be expressed as $$ \begin{align} \frac{A}{\rho}\frac{\partial p}{\partial x} = \frac{\partial}{\partial x}\left( \frac{f A^{3/2}}{3 \rho A_0^{1/2}} \right) + \left( \frac{2 A^{3/2}}{3 \rho A_0^{1/2}} - \frac{A}{\rho} \right)\frac{\partial f}{\partial x} - \frac{f A^{3/2}}{3 \rho A_0^{3/2}}\frac{\partial A_0}{\partial x} \tag{43} \end{align} $$
Thus, the system can now be expressed in conservational form as $$ \begin{aligned} \frac{\partial}{\partial t} \underbrace{\left( \begin{matrix} A \\ Q \end{matrix} \right)}_{U} + \frac{\partial}{\partial x}\underbrace{\left( \begin{matrix} Q \\ \frac{Q^2}{A} + \frac{f A^{3/2}}{3\rho A_{0}^{1/2}} \end{matrix} \right)}_{F} = \underbrace{\left( \begin{matrix} 0 \\ - \frac{2 \pi \nu q R}{\delta A} + \frac{1}{\rho}\left( \frac{2\sqrt{\pi} f A^{3/2}}{3 A_0} - \left( \frac{2 A ^{3/2}}{3A_0^{1/2}} - A \right)\frac{\partial f}{\partial r_0} \right)\frac{dr_0}{dx} \end{matrix} \right)}_{S} \end{aligned} $$
with \( U \), \( F \) and \( S \) the conservative variable, the flux and the source, respectively.
Characteristic equations
Equations in (42) form a system of hyperbolic PDEs that can be analysed with the method of characteristics. In particular, system (42) can be expressed in non-conservative form as [4] $$ \begin{align} \tag{44} \frac{\partial \mathbf{U}}{\partial t} + \underbrace{\frac{\partial \mathbf{F}}{\partial \mathbf{U}}}_{\mathbf{H}}\frac{\partial \mathbf{U}}{\partial x} = \mathbf{S} \end{align} $$
where $$ \mathbf{H} = \left[ \begin{matrix} 0 & 1 \\ -\frac{Q^2}{A^2} + c^2 & \frac{2Q}{A} \end{matrix} \right] $$
and \( c \) the wave speed $$ \begin{align} \tag{45} c = \sqrt{\frac{A}{\rho}\frac{\partial p}{\partial A}} \end{align} $$
The above system can be analysed using Riemann's method of characteristics [6]. Since \( A > 0 \) and \( \frac{A}{\rho}\frac{\partial p}{\partial A} > 0 \) in normal physiological conditions, \( H \) has two real and distinct eigenvalues (which are given by \( \left| \Lambda \mathbf{I} - \mathbf{H} \right| = 0 \)), \( \lambda_{f, b} = u \pm c \), where \( c = \sqrt{\frac{A}{\rho}\frac{\partial p}{\partial A}} \), is the pulse wave speed of the system.
Linear algebra shows \( \mathbf{H} \) to be diagonisable in the form \( \mathbf{H} = \mathbf{L}^{-1} \mathbf{\Lambda} \mathbf{L} \). The columns of \( \mathbf{L} \) are the left eigenvectors of \( \mathbf{H} \) (by solving the system \( L_i (\lambda_i \mathbf{I} - \mathbf{H} ) = 0) \)). Thus, if we left multiply system (44) by \( \mathbf{L} \), the above system becomes $$ \begin{align} \mathbf{L} \frac{\partial \mathbf{U}}{\partial t} + \mathbf{\Lambda} \mathbf{L} \frac{\partial \mathbf{w}}{\partial x} = \mathbf{LS} \tag{46} \end{align} $$
Taking \( \frac{\partial \mathbf{W}}{\partial \mathbf{U}} = \mathbf{L} \), where \( \mathbf{W}=[W_f, W_b]^T \) is the vector of characteristics (or Riemann) variables, the above system can be transformed to $$ \begin{align} \tag{47} \frac{\partial \mathbf{W}}{\partial t} + \mathbf{\Lambda} \frac{\partial \mathbf{W}}{\partial x} = \mathbf{LS} \end{align} $$
For any path \( x= \hat{x}(t) \) in the \( (x,t) \) space, the total derivative of \( \mathbf{W} \) can be expressed as $$ \begin{align} \tag{48} \frac{d\mathbf{W}(\hat{x}(t), t)}{dt} = \frac{\partial \mathbf{W}}{\partial t} + \frac{d\hat{x}}{dt} \mathbf{I} \frac{\partial \mathbf{W}}{\partial \hat{x}} \end{align} $$
From Eqs. (47) and (48), the total derivative of \( \mathbf{W} \) can be expressed as $$ \begin{align} \frac{d\mathbf{W}(\hat{x}(t), t)}{dt} = \mathbf{LS} \tag{49} \end{align} $$
The characteristic variables \( W_f \) and \( W_b \) can be determined by integrating \( \frac{\partial \mathbf{W}}{\partial U} = L \) component wise (see Sherwin et al. [7]) and they are expressed as $$ \begin{align} \tag{50} W_{f, b} = \int \frac{\partial \mathbf{W}}{\partial \mathbf{U}} = \int du \pm \int\frac{c}{A}dA = \frac{Q}{A} \pm 4c \end{align} $$
Therefore, for any point \( (X, T) \) in \( (x, t) \) space, there are two characteristic paths, \( C_f \) and \( C_b \), defined as $$ \begin{align} C_{f, b} \equiv \frac{d \hat{x}_{f, b}}{dt} = \lambda_{f, b} = U \pm c \tag{51} \end{align} $$
Along these paths, Riemann variables, \( W_f \) and \( W_b \), propagate at speeds \( \lambda_f \) and \( \lambda_b \), respectively. Under physiologic conditions, \( \lambda_f > 0 \) and \( \lambda_b < 0 \) (the flow is subcritical). Thus, the wave speed, \( c \), is generally much higher than the maximum convective velocity \( U \).
Finally, \( q \) and \( A \) variables can be calculated from the Riemann variables (by adding or subtracting the two equations in (50)) as $$ \begin{align} \tag{52} A^c = A_0 \left( \frac{W_f - W_b}{4} \right)^4 \left( \frac{2 \rho}{f} \right) \end{align} $$
and $$ \begin{align} \tag{53} Q^c = A^c \left( \frac{W_f + W_b}{2} \right) \end{align} $$
The above characteristic analysis can be interpreted as that pressure and velocity \textit{wavefronts}\footnote{\textit{wavefront} refers to an infinitesimal change in one of the Riemann variables \( W_f \) or \( W_b \) [8]} propagate forwards and backwards at a speed of \( u + c \) and \( u - c \), respectively. Forward running wavefronts are initially generated by the contraction of the heart and then they propagate to the peripheral arteries from the aorta. At braching points and/or in topologies with different mechanical properties (e.g. stenosis, aneurysms, etc), wavefronts experience reflection; as a result, they travel back to the heart. However, these reflections as they travel back might be re-reflected due to discontinuities and parts of these wavefronts are then travel towards to the periphery. This process is repeated continuously during the heart cycle and is more evident to more peripheral sites which the so-called "wave-trapping" is experienced (see wave-trapping: http://www.bg.ic.ac.uk/research/k.parker/wave_intensity_web/wia-6-2.html). Therefore, pressure and flow measured at any point in the arterial tree can be explained as the combination of many forward and backward running wavefronts. More detail on this can be found in the wave decomposition section.
Wave decomposition
According to Eq (50), the variation of the Riemann variables \( dW_{f, b} \) is related to the variation of area and velocity. From Eq. (45), the wave speed can also be expressed as \( c = (A/\rho c)(dp/dA) \). Thus, Eq. (50) can produce $$ \begin{align} \tag{54} dW_{f, b} = du \pm \frac{dp}{\rho c} \end{align} $$
\( dp \) and \( du \) have contributions from forward and backward components. Therefore, these terms can be expressed as $$ \begin{align} dp = dp_f + dp_b, \quad du = du_f + du_b \tag{55} \end{align} $$
which in turn, the forward and backward variation of can be calculated as $$ \begin{align} dp_{f, b} = \frac{1}{2}\left( dp \pm \rho c du \right), \quad du_{f, b} = \pm \frac{1}{2}\left( \frac{dp}{\rho c} \pm du \right) \tag{56} \end{align} $$
Finally, the forward and backward components of pressure an velocity waveforms can be determined by adding the instantaneous differences of the above equation over the time period \( t \), as $$ \begin{align} p_f(t) = p_0 + \sum_0^t dp_f, \quad p_b(t) = \sum_0^t dp_b \tag{57} \end{align} $$ $$ \begin{align} u_f(t) = \sum_0^t du_f, \quad u_b(t) = \sum_0^t du_b \tag{58} \end{align} $$
where \( p_0 \) is the integration constant for the pressure. The integration constant for velocity is considered to be zero. In Fig. (1), a decomposed pulse wave proximal to brachial artery is depicted.
Figure 1: Wave decomposition in proximal brachial artery.
Intensity analysis
Solving eq. (54) for \( dp(t) \) and \( du(t) \) yields $$ \begin{align} dp = \frac{\rho c}{2}\left( dW_f - dW_b \right), \quad du = \frac{1}{2}\left( dW_f + dW_b \right) \tag{59} \end{align} $$
Wave intensity is defined as $$ \begin{align} dI = dp du = \frac{\rho c}{4}\left( (dW_f)^2 - (dW_b)^2 \right) \tag{60} \end{align} $$
which is the flux of energy per unit area carried by wavefronts as they propagate through a point x. Wave intensity measures the importance of \( p \) and \( u \) in forward and backward directions. An example of the intensity waveform along with its decomposition into forward and backward components are depicted in Fig. (2).
| \( dI \) | \( dp \) | \( du \) | type of wave |
| > 0 | > 0 | > 0 | FCW |
| > 0 | < 0 | < 0 | FEW |
| < 0 | > 0 | < 0 | BCW |
| < 0 | < 0 | > 0 | BEW |
Figure 2: Intensity analysis in proximal brachial artery.
