Acta Polytechnica doi:10.14311/AP.2016.56.0099 Acta Polytechnica 56(2):99–105, 2016 © Czech Technical University in Prague, 2016 available online at http://ojs.cvut.cz/ojs/index.php/ap INTEGRAL METHODS FOR DESCRIBING PULSATILE FLOW David Hromádkaa,∗, Hynek Chlupa, Rudolf Žitnýb a Department of Mechanics Biomechanics and Mechatronics, Czech Technical University in Prague, Technická 4, 166 07 Praha 6, Czech Republic b Department of Process Engineering, Czech Technical University in Prague Technická 4, 166 07 Praha 6, Czech Republic ∗ corresponding author: david.hromadka@fs.cvut.cz Abstract. This paper presents an approximate solution of the pulsatile flow of a Newtonian fluid in the laminar flow regime in a rigid tube of constant diameter. The model is represented by two ordinary differential equations. The first equation describes the time evolution of the total flow rate, and the second equation characterizes the reverse flow near the wall. These equations are derived from the momentum balance equation and from the kinetic energy equation, respectively. The accuracy of the derived equations is compared with a solution in which the finite difference method is applied to a partial differential equation. Keywords: pulsatile flow, integral energy, integral momentum, rigid pipe. 1. Introduction Pulsatile flow is a time-dependent flow that consists of a constant part (Hagen-Poiseuille flow, with a parabolic velocity profile, fully controlled by the vis- cous forces) and an oscillatory part, which is controlled by the viscous forces at the wall and by the inertial forces at the core [1]. Pilot studies of time-dependent laminar flow were reported by [2], and were followed by [3] and [4]. Sexl [2] developed an analytical solution to the momentum equation of Newtonian fluid flow subjected to a pressure gradient dp/dx = Ceiωt. Szy- manski [3] employed a similar technique, and obtained an analytical solution to the momentum equation for the flow driven by a pressure gradient dp/dx = 0 for t ≤ 0 and dp/dx = C for t > 0, where C was con- stant. An extended review of the literature dealing with laminar pulsatile flow, which presents the current state of theoretical and experimental knowledge was thoroughly drafted by [5] and supplemented in [6]. The most cited work on linearization of the basic equations [7] was carried out by Womersley in the late 1950’s, and was published in a series of papers [4], [8] and [9], [10]. Most of these materials appear in a comprehensive report [11]. Another option for solving non-stationary flow is to employ the integral momentum and integral energy method, which leads to a different form of equations. In this method, it is not the Navier-Stokes equations themselves that are solved, but the integrated form of the stream wise momentum or energy equation. The form of the velocity profiles is assumed a priori with unknown coefficients characterizing the profile [12]. A list of integral methods that are used, variants of velocity profiles, the ability to capture the reverse flow and specific features of the method are presented in (Tab. 1). The methods in the list describe unidirec- tional the pulsatile flow and the impulsively-initiated flow of a Newtonian fluid within a rigid tube of con- stant diameter, unless otherwise stated in the column “Specific”. To the best of the author’s knowledge, there is no in- tegral energy solution for pulsatile flow that takes into account the reverse flow. Only some of the listed meth- ods are able to solve both, the periodic flow and the impulsively-initiated flows for arbitrary initial condi- tions (using time marching algorithms). Additionally, there is a controversy over the comparative accuracy of the energy method and momentum integral meth- ods; the results depend on the specific test and on the selected accuracy criterion. For example, Elkouh [21] tested his method in the special case of a sinusoidally oscillating Newtonian fluid (the method did not take into account the reverse flow near the wall). A com- parison with the exact solution [22] showed, that the integral energy method is more accurate than the in- tegral momentum method. The methods described in this paper anticipate quite opposite conclusions. The method proposed in our paper was applied in the newly developed experimental method for identify- ing viscoelastic properties of blood vessels and grafts using a transient water hammer experiment [24]. The experimental setup consists of the tested elastic tubu- lar sample connected to a long glass capillary filled with an oscillating column of water. When the stan- dard balance momentum equation with a parabolic ve- locity profile [25], [26] is used to describe the pulsation of the column of water, there is an error in the pre- dicted frequency. The error is reduced by using a new model, which replaces the standard balance momen- tum equation by two ordinary differential equations: one for the total flow rate and one for the reverse flow rate. The model is derived from the momentum bal- ance and from the macroscopic kinetic energy balance. 99 http://dx.doi.org/10.14311/AP.2016.56.0099 http://ojs.cvut.cz/ojs/index.php/ap D. Hromádka, H. Chlup, R. Žitný Acta Polytechnica Author Solution method Velocity profile form Reverseflow Specific [12] Integral momentum u(r,t) = U(t) 4∑ i=0 ai ri Ri for a1 = 0 yes Flow in a tube with and without longitudinal vibration of the wall compared with [11] exact solution for α = 0–10 [13] Integral momentum u(r,t) = U(t) 4∑ i=0 ai ri Ri for a1 = a3 = 0 yes Flow of Non-Newtonian fluid with vibrating wall [14] Integral momentum u(r,t) = U(t) 4∑ i=0 ai ri Ri for a1 = 0 yes Flow of micropolar fluid with and without longitudinal vibration of the wall [15] Combination of integral energy and integral momentum u(r) = U 4∑ i=0 ai ri Ri for a1 = a3 = 0 and flat profile yes Steady flow through rigid stenosis, results were compared with an experiment [16] Integral momentum u(r) = U 4∑ i=0 ai ri Ri for a0 = 0 yes Steady flow through rigid stenosis, results were compared with experiment [17] [18] Integral momentum u∗(r∗,x∗, t∗) = 3∑ i=0 Ui(x∗, t∗)r∗i for U0 = U1 = 0 yes Flow through compliant and rigid tube. Solution in rigid tube was compared with [11] for α = 1, 3, 4.72, 6.67, 15 [19] Integral momentum u(r,t) = eiωta0 ( 1− r 2 R2 )( a1 − r 2 R2 ) yes Solution only for harmonic P, compared with [20] exact solution for α= 2, 4, 6 [21] Integral momentum,integral energy u(r,t) = 2Ū(t) 3n+1 n+1 ( 1 − r n+1 n R n+1 n ) tested only for the special case n = 1 no Flow of power law and Newtonian fluid compared with [22] exact solution for α = 0–10 [23] Integral momentum u(r,t) = 2U(t) ( 1 − r 2 R2 ) no Method proposed compared with [22] exact solution for α = 0–9 Table 1. List of integral methods used for describing pulsatile flow. 2. Methods 2.1. Mathematical model – Weak formulation The exact solution of pulsatile flow in a rigid tube of constant radius R is fully described by the momentum balance in the x axial direction, assuming rotational symmetry and spatially fully-developed flow of a non- compressible fluid ρ ∂u ∂t = − ∂p ∂x + µ 1 r ∂ ∂r ( r ∂u ∂r ) , (1) where ρ is constant density and µ is dynamic viscosity. Pressure p is independent of the radial coordinate, as follows from the momentum balance in the radial direction. The pressure gradient can be considered as a known (and arbitrary) function P(t) P(t) = ∂p ∂x . (2) The equation (1) multiplied by velocity u represents the kinetic energy balance, which forms the basis of the mechanical energy equation ρ 1 2 ∂u2 ∂t = −uP + uµ 1 r ∂ ∂r ( r ∂u ∂r ) . (3) An exact solution of (1), (3) has to provide the same results, but the results of an approximate solution may differ slightly. An approximation of the velocity profile u(r,t) can be represented as a linear combination of polynomial basis functions u(r,t) = n∑ i=1 Ui(t)Ni(r) = n∑ i=1 Ui(t) ( r2(i−1) R2(i−1) − r2i R2i ) . (4) The term Ui(t) is the amplitude of the velocity profile, and the symmetric basis functions Ni(r) satisfy the boundary conditions (u(R,t) = 0,∂u(0, t)/∂r = 0). 100 vol. 56 no. 2/2016 Integral Methods for Describing Pulsatile Flow An approximation of the velocity profile for n = 1 corresponds to a parabolic function (Hagen Poiseuille) and further terms with a higher power exponent can be called corrective functions that provide a description of more extrema (for example the second basis function N2(r) represents a typical radial velocity profile with backflow at the wall). First, we derive approximate solutions (5), (6) from the balance momentum equation (1) and from the kinetic energy equation (3) respectively. The Galerkin weighted residual method is applied. The integration of the residual res and the weight function wj is over area 2πr dr, and for the momentum balance results in R∫ 0 res · wj dr = R∫ 0 [ ρ ∂u ∂t + P − µ 1 r ∂ ∂r ( r ∂u ∂r )] Nj (r)2πr dr = ρ n∑ i=1 ∂Ui (t) ∂t R∫ 0 Ni(r)Nj (r)2πr dr ︸ ︷︷ ︸ M b − µ n∑ i=1 Ui(t) R∫ 0 Nj (r) 1 r ∂ ∂r ( r ∂Ni(r) ∂r ) 2πr dr ︸ ︷︷ ︸ Kb + P R∫ 0 Nj (r)2πr dr ︸ ︷︷ ︸ nb = 0, (5) R∫ 0 [ ρ 1 2 ∂u2 ∂t + uP − uµ 1 r ∂ ∂r ( r ∂u ∂r )] Nj (r)2πr dr = ρ n∑ i=1 ∂U2i (t) ∂t R∫ 0 1 2 N2i (r)Nj (r)2πr dr ︸ ︷︷ ︸ M k − µ n∑ i=1 U2i (t) R∫ 0 Ni(r)Nj (r) 1 r ∂ ∂r ( r ∂Ni(r) ∂r ) 2πr dr ︸ ︷︷ ︸ Kk + n∑ i=1 Ui(t)P R∫ 0 Ni(r)Nj (r)2πr dr ︸ ︷︷ ︸ N k = 0. (6) Subscript i is the summation index, and index j cor- responds to the j-th weight function (corresponding to the j-th equation). It is more convenient to work with flow rates. We therefore define the total flow rate within the tube q(t) = 2π n∑ i=1 Ui(t) R∫ 0 rNi(r) dr = πR2 n∑ i=1 Ui(t) i(i + 1) = q1(t)+q2(t) + · · · + qn(t).res. (7) Subscript i = 1 corresponds to the Poiseuille flow, while indices i ≥ 2 to n denote correction functions that are able to describe more extrema of the velocity profile qi(t) = 2πUi(t) R∫ 0 rNi(r) dr = πR2 Ui(t) i(i + 1) . (8) It is now simple to express the velocity Ui (t) as a function of the flow rate  U1 (t) Ui (t) ... Un (t)   ︸ ︷︷ ︸ u = 1 πR2   2 −2 −2 −2 0 i(i+1) · · · 0 ... 0 ... 0 0 0 0 n(n+1)   ︸ ︷︷ ︸ A   q (t) qi (t) ... qn (t)   ︸ ︷︷ ︸ q . (9) Equations (5) and (6) can be rewritten in terms of the flow rate in matrix notation as ρ dq dt + M −1b A −1nbP −µM −1b A −1KbAq = 0, (10) ρ d dt (q2)+M −1k A −2N kAqP −µA−2M −1k KkA 2q2 = 0. (11) Let us restrict the approximation to the first term in the polynomial series (4) (parabolic velocity profile). The momentum balance reduces into the equation ρ dq dt + 3 4 PR2π + 24 4 µ R2 q = 0 (12) and the mechanical energy equation formulated from (3) ρ dq dt + 2 3 PR2π + 16 3 µ R2 q = 0. (13) These two equations describe the same physical prob- lem. However, they are not identical as they differ from each other by the coefficients. The equations are identical only for stationary flow, and both lead to Hagen-Poiseuille flow. The reader has to ask himself which of the equations – (12) or (13) – provides a better description of the real movement of a fluid. However, the answer to this question is left until an analysis of the test examples has been made. We consider the velocity profile as the fourth order poly- nomial function (the first basis function describes the Poiseuille flow, and the second basis function describes 101 D. Hromádka, H. Chlup, R. Žitný Acta Polytechnica the reverse flow, the backflow in the vicinity of the wall) u(r,t) = U1(t) ( 1 − r2 R2 ) + U2(t) ( r2 R2 − r4 R4 ) . (14) Equations (9), (10) were applied to (14), and this leads to a system of ODE’s. After a simple manipulation we obtain a standard formulation of the momentum balance for the flow rates ρ   dq dt dq2 dt   + πR2   8 9 5 9  P + µ R2   64 9 80 9 40 9 320 9  (q q2 ) = 0. (15) And by applying (9), (11), respectively, to the velocity profile mentioned above, we obtain the formulation of kinetic energy for the flow rates ρ   dq dt dq2 dt   + πR2   1 6 8q2q + 12q22 − 35q2 2q22 + 2q2q − 7q2 7 12 −5q2 + 2q22 2q22 + 2q2q − 7q2  P+ µ R2   4 3 −20q2q2 +24q32 +4q22q−35q3 2q22 +2q2q−7q2 14 3 −5q3 +2q22q−40q2q2 +16q32 2q22 +2q2q−7q2   = 0. (16) The reader will notice that the system of equations (15) and (16) deals with greater differences than for (12) and (13). Ordinary differential equations (15) are linear, while (16) is a system of nonlinear differ- ential equations. Systems (12), (13), or (15), (16) can be solved with the Matlab ODE toolbox for an arbitrary time course of P(t) (continuous as well as discontinuous pressure gradients) and for arbitrary initial conditions prescribed in the form of the total flow rate q and the reverse flow rate q2 at time zero. 2.2. Test example Differences between the solution of the flow rate re- alized by the finite difference method qfd with a very fine mesh (an implicit method with a central differ- ence scheme), the standard balance momentum (12), the mechanical energy equation (13) and the pro- posed two-equation approximations (15), (16) were tested using the model with harmonic driving force (pressure gradient P(t)), consisting of a stationary part pstat = −100 Pa m−1 and an oscillatory part pamp = −500 Pa m−1 in a rigid tube of constant ra- dius R = 2 mm P(t) = pstat + pamp sin(ωt). (17) Stationary flow is assumed at time t = 0 as initial conditions, therefore q(0) = − πR4 8µ P(0), q2 (0) = 0. (18) The frequency ω can be expressed in a dimensionless form either Strouhal or Womersley number α α = R √ ωρ µ . (19) 2.3. Determining the deviation The deviation (20) is computed as the standard Eu- clidean norm of two flow functions, which is calculated from the finite difference method and from the flow rate using the equation (12), (15), (13), (16), respec- tively, and which are integrated over time with the period T = 2π/ω Error =   1 T T∫ 0 (q(t) − q(t)fd) 2 max (q(t)2fd) dt   1 2 . (20) The error defined in this way is dimensionless, and it is a suitable measure for deviations in the total calculated flow rates and also in the back flow rates. 2.4. Wall shear stress The wall shear stress is calculated from the following equation τ = µ ∂u ∂r . (21) 3. Results The deviation of the suggested balance momentum (12),(15), the energy equations (13), (16) is demon- strated by (Fig. 1), which illustrates the deviation from the exact solution computed by the equation (20). 0 0.01 0.02 0.03 0.04 0.05 0.06 0 2 4 6 8 10 12 14 16 18 20 E r r o r [− ] α[−] Figure 1. The graph presents the deviation of flow rate q from qfd. The black curves correspond to the formulation from the kinetic energy balance (13), (16) while the grey curves correlate with the weak solution of the balance momentum (12), (15). The solid curves are consistent with the parabolic velocity profile ap- proximation and the dashed curves are equivalent to the approximation of the velocity profile with a fourth order polynomial function. Three representative samples of the evolution of the velocity profile and the course of the wall shear 102 vol. 56 no. 2/2016 Integral Methods for Describing Pulsatile Flow stress are shown in the figures bellow. Pairs of veloc- ity profiles are presented in (Fig. 2). A Womersley number is assigned to each pair of velocity profiles in two different time points t1, t2 and also for the wall shear stress. The results of the wall shear stress can be found in (Fig. 3) below. -0.15 -0.1 -0.05 0 0.05 -2-1.5-1-0.5 0 0.5 1 1.5 2 u [m s − 1 ] r[mm] α = 3 t1 = 0.7T t2 = 0.9T 0 0.05 0.1 0.15 -2-1.5-1-0.5 0 0.5 1 1.5 2 u [m s − 1 ] r[mm] α = 5 t1 = 0.7T t2 = 0.9T 0 0.05 0.1 0.15 -2-1.5-1-0.5 0 0.5 1 1.5 2 u [m s − 1 ] r[mm] α = 7 t1 = 0.7T t2 = 0.9T Figure 2. The black lines are representatives of the velocity profiles obtained from the finite difference method on the right side of the figures. The dashed lines were obtained by calculating (16) (black color), (15) (grey color). 4. Conclusion The method proposed here is suitable especially for the situation where the pressure gradient cannot be decomposed into Fourier series, which corresponds to [24]. The exact solution through the Bessel function is therefore not applicable. The effect of assuming a parabolic profile and applying the Galerkin method is shown in (Fig. 1). The maximum error when us- ing the momentum integral method is 2.8 percent. The maximum error when using the energy integral method is 5.5 percent. The discrepancy between Elk- ouh’s results, which assume the steady state of the velocity profile, and our results is caused by the use a different method for finding the solution (in our case, the Galerkin method), by the specific test, and by the selected accuracy criterion. Elkouh reported that better results were obtained from the integral energy method with the well-known steady state velocity profile (without reverse flow) for α = 0–10 [21]. The method proposed here is -2 0 2 4 0 0.5 1 1.5 2 2.5 τ (t ) τ (0 )[ − ] t[s] α = 3 0 1 2 3 0 0.2 0.4 0.6 0.8 1 τ (t ) τ (0 )[ − ] t[s] α = 5 0 1 2 0 0.1 0.2 0.3 0.4 0.5 τ (t ) τ (0 )[ − ] t[s] α = 7 Figure 3. The black lines correspond to the shear stress at the wall calculated using the finite difference method. The dashed lines mark the wall shear stress of the approximate solution. The color designation for the approximate solution is the same as in previous figures. designed to take the reverse flow into account, and it investigates the error of the integral energy method and the integral momentum method. The maximum error over a cycle was 0.3 percent for the integral momentum solution while the maximum error caused using the integral energy method was 1.3 percent. The higher error of the system of mechanical energy equations derived from the integral energy (16) is caused by nonlinearities. The integral solution (15) is in excellent agreement with the finite difference solution in the region α = 1–20 which is a part of the intermediate region of pulsatile flow [5] , [6]. Another type of gradient P (a symmetric triangular wave) was applied to the system of equations (15), (16) and also to the finite difference method. The same amplitude of the stationary and nonstationary part of the pressure gradient was used as in the test example (17). Similar error was calculated as in (Fig. 1) based on (20). The maximum error from the flow calculation (15), (16) in comparison with the finite difference solution, was 1.4, 2.4 percent, respectively. The evolution of the approximate velocity profiles is in good agreement with the velocity profile calculated from the finite difference method (Fig. 2). The wall shear stress calculated from the momen- tum balance (15) describes the shear stress at the 103 D. Hromádka, H. Chlup, R. Žitný Acta Polytechnica wall more precisely than the shear stress that was evaluated from the mechanical energy equation (16). Other variational methods were tested (least square integral energy and momentum, an expert estimate of the weight function with integration of momentum, and also the energy equation) together with various integrations of differential equations (over the radius and over the area). The lowest error found was using the Galerkin integral momentum method with inte- gration over the area of rigid pipe. This method has therefore been presented in this paper. 5. Appendix For illustration, let us derive the balance momentum with the backflow (15). We recall equations (1), (5), (14) and we obtain a set of integral equations (22), (23) R∫ 0 { ρ [ dU1(t) dt ( 1 − r2 R2 ) + dU2(t) dt ( r2 R2 − r4 R4 )] + P − µ r [ −U1(t) 2r R2 + U2(t) ( 2r R2 − 4r3 R3 )] − µ [ −U1(t) 2 R2 + U2(t) ( 2 R2 − 12r2 R4 )]} ( 1 − r2 R2 ) 2πr dr = 0, (22) R∫ 0 { ρ [ dU1(t) dt ( 1 − r2 R2 ) + dU2(t) dt ( r2 R2 − r4 R4 )] + P − µ r [ −U1(t) 2r R2 + U2(t) ( 2r R2 − 4r3 R3 )] − µ [ −U1(t) 2 R2 + U2(t) ( 2 R2 − 12r2 R4 )]} ( r2 R2 − r4 R4 ) 2πr dr = 0. (23) The inner product of differential equations and weighted basis functions integrated over the area yield two ordinary differential equations (24), (25) for un- known center-line velocity U1(t), U2(t) 1 3 R2ρ dU1(t) dt + 1 12 R2ρ dU2(t) dt + 2µU1 (t) + 2 3 µU2(t) + 1 2 R2P = 0, (24) 1 12 R2ρ dU1(t) dt + 1 30 R2ρ dU2(t) dt + 2 3 µU1(t) + 2 3 µU2(t) + 1 6 R2P = 0. (25) We can formulate the center-line velocity in terms of flowrates. This can be done easily through (9) U1(t) = 1 πR2 ( 2q(t) − 2q2(t) ) , (26) U2(t) = 1 πR2 ( 6q2(t) ) . (27) Combining (24),(25) with (26), (27) leads to a set of a differential equations for describing the non-stationary flowrates within a rigid pipe of constant diameter ρR2 ( 4 −1 5 1 ) dq dt dq2 dt  +πR4(35 ) P +µ ( 24 0 40 80 )( q q2 ) = 0. (28) Simple manipulation of (28) yields the standard set of equations (15). Equations (16) are derived in a similar manner. List of symbols α Womersley number [–] ρ Water density [kg m−3] τ Wall shear stress [Pa] µ Dynamic viscosity [Pa s] ω Angular frequency [rad s−1] i i-th index [–] i Imaginary unit [–] n n-th index [–] Ni Basis function [–] p Pressure [Pa] P Pressure gradient [Pa m−1] pamp Amplitude of the nonstationary part of gradient pressure [Pa m−1] q Total flowrate [m3 s−1] qi Corrective flowrate [m3 s−1] qfd Total flowrate calculated from finite difference method [m3 s−1] r Radial coordinate [m] r∗ Dimensionless radial coordinate [–] R Radius of the tube [m] t Time [s] t∗ Dimensionless time [–] u Fluid velocity [m s−1] u∗ Dimensionless velocity [–] Ui Velocity amplitude [m s−1] Ū Average velocity amplitude [m s−1] wj j-th weight function [m] x Axial coordinate [m] x∗ Dimensionless axial coordinate [–] res Residual [Pa m] Error Deviation of the flowrate [–] A Matrix flowrate coefficients [–] Kb Stiffnes matrix of the momentum balance [–] Kk Stiffnes matrix of the kinetic balance [–] M b Mass matrix of the momentum balance [m2] M k Mass matrix of the kinetic balance [m2] N b Coefficients of driving force [m2] N k Coefficients of driving force [m2] q Flowrate vector [m3 s−1] 104 vol. 56 no. 2/2016 Integral Methods for Describing Pulsatile Flow Acknowledgements This work has been supported by Czech Technical Univer- sity in Prague grant No. SGS13/176/OHK2/3T/12 and by Ministry of Health project Nos. CR NT 13302 and 15-27941A. References [1] M. Zamir. The Physics of Pulsatile Flow. Springer Science Business Media, 2000. doi:10.1007/978-1-4612-1282-9. [2] T. Sexl. Über den von e. g. richardson entdeckten annulareffekt. Zeitschrift für Physik 61(5-6):349–362, 1930. doi:10.1007/bf01340631. [3] P. J. Szymanski. Quelques solutions exactes des équations de l’hydrodynamique du fluide visqueux dans le cas d’un tube cylindrique. J Math Pures Appl 9(11):67–107, 1932. [4] J. R. Womersley. Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known. The Journal of Physiology 127(3):553– 563, 1955. doi:10.1113/jphysiol.1955.sp005276. [5] M. Y. Gundogdu, M. O. Carpinlioglu. Present state of art on pulsatile flow theory. Part 1. Laminar and transitional flow regimes. JSME International Journal Series B 42(3):384–397, 1999. doi:10.1299/jsmeb.42.384. [6] M. O. Carpinlioglu, M. Y. Gundogdu. A critical review on pulsatile pipe flow studies directing towards future research topics. Flow Measurement and Instrumentation 12(3):163–174, 2001. doi:10.1016/s0955-5986(01)00020-6. [7] G. Rudinger. Review of current mathematical methods for the analysis of blood flow. ASME Biomedical Fluid Mechanics Symposium pp. 1–33, 1966. [8] J. R. Womersley. Oscillatory flow in arteries: the constrained elastic tube as a model of arterial flow and pulse transmission. Physics in Medicine and Biology 2(2):178–187, 1957. doi:10.1088/0031-9155/2/2/305. [9] J. R. Womersley. Oscillatory flow in arteries. III: Flow and pulse-velocity formulae for a liquid whose viscosity varies with frequency. Physics in Medicine and Biology 2(4):374–382, 1958. doi:10.1088/0031-9155/2/4/307. [10] J. R. Womersley. Oscillatory flow in arteries. II: The reflection of the pulse wave at junctions and rigid inserts in the arterial system. Physics in Medicine and Biology 2(4):313–323, 1958. doi:10.1088/0031-9155/2/4/301. [11] J. R. Womersley. An elastic tube theory of pulse transmission and oscillatory flow in mammalian arteries. Tech. Rep. TR56–614, WADG Technical Report, 1957. [12] L. E. Hooks, R. M. Nerem, T. J. Benson. A momentum integral solution for pulsatile flow in a rigid tube with and without longitudinal vibration. International Journal of Engineering Science 10(12):989– 1007, 1972. doi:10.1016/0020-7225(72)90021-3. [13] J. Misra, B. Kar. Unsteady flow of blood through arteries in vibration environments. Mathematical and Computer Modelling 13(4):7–17, 1990. doi:10.1016/0895-7177(90)90049-s. [14] S. Parvathamma, R. Devanathan. Microcontinuum approach to the pulsatile flow in tubes with and without longitudinal vibration. Bulletin of Mathematical Biology 45(5):721–737, 1983. doi:10.1007/bf02460045. [15] B. E. Morgan, D. F. Young. An intergral method for the analysis of flow in arterial stenoses. Bulletin of Mathematical Biology 36(1):39–53, 1974. doi:10.1007/bf02461189. [16] J. H. Forrester, D. F. Young. Flow through a converging-diverging tube and its implications in occlusive vascular disease — I. Journal of Biomechanics 3(3):297–305, 1970. doi:10.1016/0021-9290(70)90031-x. [17] J. H. Forrester, D. F. Young. Flow through a converging-diverging tube and its implications in occlusive vascular disease — II. Journal of Biomechanics 3(3):307–316, 1970. doi:10.1016/0021-9290(70)90032-1. [18] F. K. Tsou, P. C. Chou, S. N. Frankel, A. W. Hahn. An integral method for the analysis of blood flow. The Bulletin of Mathematical Biophysics 33(1):117–128, 1971. doi:10.1007/bf02476669. [19] E. F. Blick, P. D. Stein. Variational solution for pulsatile flow in rigid tubes. Journal of Biomechanical Engineering 106(1):89, 1984. doi:10.1115/1.3138464. [20] P. Lambossy. Oscillations forcées d’un liquide imcompressible et visqueux dans un tube rigide et horizontal. Calcul de la force de trottement. Helv Physica Acta 25, 1957. [21] A. F. Elkouh. Approximate solution for pulsatile laminar flow in a circular rigid tube. Journal of Fluids Engineering 100(1):131, 1978. doi:10.1115/1.3448588. [22] S. Uchida. The pulsating viscous flow superposed on the steady laminar motion of incompressible fluid in a circular pipe. Journal of Applied Mathematics and Physics (ZAMP) 7(5):403–422, 1956. doi:10.1007/bf01606327. [23] J. McFeeley, R. Patel, K. Jolls. Approximate low-frequency solution for pulsatile laminar flow in a tube. Chemical Engineering Science 28(11):2105–2107, 1973. doi:10.1016/0009-2509(73)85058-4. [24] D. Hromádka, H. Chlup, R. Žitný. Identification of relaxation parameter of a physical model of vein from fluid transient experiment. EPJ Web of Conferences 67:02039, 2014. doi:10.1051/epjconf/20146702039. [25] R. Bird, W. Stewart, E. Lightfoot. Transport Phenomena. Wiley International edition. Wiley, 2007. [26] I. Idelchik. Handbook of Hydraulic Resistance. Jaico Publishing House, 2005. 105 http://dx.doi.org/10.1007/978-1-4612-1282-9 http://dx.doi.org/10.1007/bf01340631 http://dx.doi.org/10.1113/jphysiol.1955.sp005276 http://dx.doi.org/10.1299/jsmeb.42.384 http://dx.doi.org/10.1016/s0955-5986(01)00020-6 http://dx.doi.org/10.1088/0031-9155/2/2/305 http://dx.doi.org/10.1088/0031-9155/2/4/307 http://dx.doi.org/10.1088/0031-9155/2/4/301 http://dx.doi.org/10.1016/0020-7225(72)90021-3 http://dx.doi.org/10.1016/0895-7177(90)90049-s http://dx.doi.org/10.1007/bf02460045 http://dx.doi.org/10.1007/bf02461189 http://dx.doi.org/10.1016/0021-9290(70)90031-x http://dx.doi.org/10.1016/0021-9290(70)90032-1 http://dx.doi.org/10.1007/bf02476669 http://dx.doi.org/10.1115/1.3138464 http://dx.doi.org/10.1115/1.3448588 http://dx.doi.org/10.1007/bf01606327 http://dx.doi.org/10.1016/0009-2509(73)85058-4 http://dx.doi.org/10.1051/epjconf/20146702039 Acta Polytechnica 56(2):99–105, 2016 1 Introduction 2 Methods 2.1 Mathematical model – Weak formulation 2.2 Test example 2.3 Determining the deviation 2.4 Wall shear stress 3 Results 4 Conclusion 5 Appendix List of symbols Acknowledgements References