CC BY
20
VISCOELASTIC EFFECTS OF BLOOD FLOW IN NONDEFORMABLE
CAPILLARY
S.N. Aristov, O.I. Skul'skiy
Institute of Continuous Media Mechanics, 1, Korolev Street, 614013, Perm, Russia
Abstract: A polydisperse structure of blood, the major portion of which is made up of erythrocytes accounts for its complicated rheological behavior. To a first approximation, the blood flow in a vessel is simulated by a stationary motion of viscoelastic fluid in a capillary with nondeformable walls. An exact analytical solution has been developed for a Poiseuille flow of viscoelastic fluid whose behavior is described by the generalized 6-constant Jeffreys model. The velocity profile is obtained in a parametric form where the velocity gradient at the capillary wall is used as a parameter. It has been demonstrated that at pressure gradients exceeding the obtained critical values the longitudinal velocity profiles show slight tangential discontinuities. The predicted phenomenon of hysteresis is shown to be associated with a monotonic decrease/increase in the pressure gradient over a certain range of parameters.
Key words: blood flow, non-Newtonian fluid, exact solution, tangential discontinuities, hysteresis
Introduction
It is a well established fact that blood is a concentrated polydisperse suspension of particles with complicated properties which are still not clearly understood. The blood flow in a vessel is accompanied by reversible aggregation and diffusion of erythrocytes as well as by orientational effects [1]. An immediate task in studying any complex medium is to decide which of its properties are most essential for constructing its simplified model and to identify the basic microscopic variables. This task suggests the introduction of state variables to provide a quantitative description of the microstructure s. Since the state of rest in a real medium generally reflects thermomechanical equilibrium of the material as a whole, it should be regarded as an isotropic one and any variation in this state is hindered by some recovery mechanism. Each of the recovery mechanisms is associated with the relaxation time X. Description of the flow at the microscopic level requires usage of the evolutionary equations relating the structure development to the mean flow and recovering force. In transition to a macroscopic level it is necessary to employ the procedure of averaging the representative element over the volume S = (s) . The general equations for a Newtonian suspended fluid are represented as [2] :
S = S(S, Vv; X), (1)
o = a (s, Vv; X), (2)
where S is the tensor parameter of order determining the suspension macrostructure, o is the stress tensor, Vv is the velocity gradient, and X is the relaxation time. The dot over the symbol denotes the objective time derivative. When deviation of the microstructure from equilibrium equations is small (1) and (2) can be simplified as:
S = aD - j3AS, (3)
o = aD + bXS, (4)
where D is the symmetric part of Vv, a,fi,a,b are material parameters.
In the case of a weak flow (|D\ << A), the approximate solution of the evolutionary equation (3) takes the form
t
S = JaD(r) exp(- pA(t - z))rfz. (5)
0
Substituting (5) into (4) gives the constitutive equation in the form consistent with linear viscoelasticity
t
ct = aD + bA^aD(z) exp(- f3A(t - z))dz . (6)
0
Thus, the models of viscoelastic fluids can be used as the first approximation for describing suspension flows including blood flow. In the present paper, the blood flow in a vessel is modeled by a stationary Poiseuille flow in a capillary with rigid walls.
The problems dealing with viscoelastic fluid behaviour, despite their long-standing investigation history are still crucial for rheology. One way to describe viscoelastic fluid properties is to construct differential models of relaxation type which are best represented by the classical Maxwell and Jeffreys models. Since the advent of these models a number of their generalized versions has been proposed, which can be differentiated by the type of time derivative used in the model equations.
The performance of different models are normally tested in the Couette-type flows, in which the flow kinematics is known and parameters to be determined are the effective viscosity and the first and the second differences of normal stresses. In the Poiseuille flows, the velocities are determined by solving boundary value problem. For this type of flows an exact integration of equations of motion yielding parabolic velocity profiles is possible only in a few particular cases. The analysis of some special cases of models has shown that a unidirectional flow in a channel can yield multivalued velocity profiles.
To our knowledge not every solution of the motion equation, be it exact or not, is actually realized in nature. The obtained velocity profiles must be dynamically and kinematically admissible, which means that they must fit the hydrodynamic equations and boundary conditions and at the same time lie within the examined geometrical region. Any type of motions generated in nature must be normally stable: small disturbances once initiated will decay with time. The challenge of making stability analysis for this class of problems is coupled with the absence of well-defined general nonlinear theory of stability for viscoelastic fluids, which complicates the selection of meaningful realistic solutions. However, earlier investigations [3-13] have thrown some light on the situation. It has been found that the flow of Newtonian fluid in a tube is unstable provided that the flow rate increases with decrease of the pressure gradient. This conclusion appears to be true for non-Newtonian fluids. The results of investigations by Truesdell and Noll [3], Coleman [4], Denn and Fosdic [5] show that there exist kinematically admissible flows which lead to negative values of the stress power. It is hardly probable that such motions will occur in nature. The investigation of the 4-constant Oldroyd fluid flow in a tube [6] has revealed nonuniqueness of the velocity gradient distribution in a radial direction. The calculated velocity profiles were either of smooth parabolic form or had an inflection with derivative discontinuity depending on the values of pressure drop and the ratio of the retardation to relaxation time. The numerical analysis indicates that in this case the flow curve shows hysteresis. The investigation of the flow around a body made by Altman and Denn [7] for a viscoelastic Maxwell-type fluid led authors to the conclusion that the condition El2 = 1 demarcates two different situations. In the subcritical case (El2 < 1), the equations are elliptic and have continuous solutions. In the
supercritical case (El2 > 1), the equations are hyperbolic and their solutions have drastic tangential discontinuities.
In the present work, our main focus is the generalized version of the Jeffreys model represented in terms of time derivatives Fcbc with arbitrarily chosen material constants a, b, c.
This version is a special case of the 8-constant Oldroyd model, in which the time derivatives of excess stresses and strain rates involve the same material constants. With the notation used in [8] a transition from the 8-constant Oldroyd model to the 6-constant Jeffreys model can be accomplished under the requirement that
/u1 = -aAj, /u0 = cX1,v1 = bX1, /u2 = -aX2, v2 = bA2.
In turn, the Maxwell, DeWitt, White-Metzner and 3 and 4-constant Oldroyd models should be considered as special cases of the 6-constant Jeffreys model.
The objectives of the present work are twofold: to find analytical solutions to the equation of momentum conservation equation for the 6-constant generalized Jeffreys model in terms of one-dimensional flow approximation; to select from the available set of solutions physically valid versions using the results of stability theory for viscoelastic flows.
Statement of the problem and exact solution
Let us consider an exact formulation for a Poiseuille flow of viscoelastic fluid in a slit described by the 6-constant Jeffreys model. The equation of motion written in terms of extra stress takes the form [9]
(7)
where T = T + (a-p)l is extra stress, T* is the stress deviator, p is the pressure, a is an arbitrary scalar.
The fluid is assumed incompressible
V • v = 0. (8)
The Jeffreys equation with the most general associated derivative retaining symmetry is written as T + X1 FabcT = ¡u(D + X2 Fabc D) and then expended to
T + A,
--W • T + T • W + a(T • D + D • T) + b(T : D)l + cDtr (T)
dt
D + ^2 (- W • D + D • W + 2a(D • D) + b(D : D)l
(9)
= 2^
where Vv = D + W is the velocity gradient, D = 1 (vvT + Vv) is the symmetric part of Vv,
W = - (VvT - Vv) is the asymmetric part of Vv.
Equation (9) contains 6 unknown material constants, which should be determined from rheological tests. We shall restrict our consideration to the case when the rheological behavior of materials is described by equation (9) with the constants meeting the requirement (a + b)(a + c) < 1. The special linear cases of general equation (9) such as the Oldroyd and Maxwell equations with lower and upper convected derivatives lead to a trivial solution with a parabolic velocity profile and will be not discussed here.
In the case of a unidirectional (vz (r), vr = v^ = 0) steady-state flow with constant coefficient of dynamic viscosity n, all components of the stress tensor are independent of z,
dp
all derivatives with respect to z, except for — = R = const, are equal to zero and the equation
dz
of motion has the form
T = 11 ^ ) ■ (10)
oz r or
after integration we have
2x w = -Rr. (11)
The components of the stress tensors are related by
т rr + т rz ^i
f ^1(1 + a + b)=|A 2 {^2 (1 + a + b), (12)
v dr ) ^ dr )
dv, V „ f dv,x2
тфф + тrz|b = |A2| ~r I b, (13)
dr ) ^ dr _
2
т zz + |(-1 + a + b) = 21 ^ I (-1 + a + b), (14)
dr ) I dr
f dvz
тrz + 2 ^ f Ifr )[(тzz - тrr ) + (Tzz + Trr )(a + С)] = I ' (15)
Making use of (12) - (14) and substituting (15) into (11) we find
dv2
(16)
2т« = 2M
fк Л1 + ^2[1 -(a + b)(a + c)]{^f
v dr ) 1 + ^2[1 -(a + b)(a + c)]fdV^2
= - Rr.
- dr )
From (16) we can readily obtain the expressions for the effective viscosity, and the first and second differences of the normal stresses
1 + 2 [1 - (a + b )(a + c)] | dVz
Л
f dvz Л I dr
= M-, (17)
V
dr ) . .9 r. ✓ 4ifdv
1 + ^21 [1 -(a + b )(a + c )]\^
dr
Ш = (тzz -тrr ) =_2M 2 )__(18)
Ш = ГИ \2 = TT72, (18)
fdfr) 1 + "12 [1 -(a + b )(a + c)]{f)
Ш =(тrr -тфф)= -M(1 +a)(^1 2) (19)
Ш = iA \2 = iA ^ . (19)
f ^ I 1 + ^2[1 -(a + b)(a + c)]| dvz dr ) I dr
In order to find solution for a viscoelastic flow in a capillary under the pressure gradient we must consider equation of motion (16) with the assumption that (a + b)(a + c) ^ 1. Otherwise we have the trivial solution with parabolic velocity profile.
r k 2
After introducing the dimensionless parameters £ = —, s = —,
rn ki
A = 71 -{a + b)(a + c), w = VzklA
r0
^k Ar
G = —1—0, we obtain from (16)
I
and new dimensionless variables q = k —^A and
dr
2q ^ = -G£. 1 + q2
(20)
Exactly the same equation (up to notations) was first derived for the 4-constant Oldroyd model in [6]. For a given dimensionless pressure drop G it describes the variation of the dimensionless velocity gradient q in a radial direction as a function of £ . The key feature of equation (20) is that due to a nonlinear character its solution is not unique as demonstrated in Fig. 1.
In the limited range of s< 1/9, the function q = ^(G£) has a region defined by a minimum and maximum of G£ where the dimensionless velocity gradient has three values. The smaller is s the wider is the region of nonuniqueness and the larger is the discontinuous jumps of the velocity gradient. This means that velocity profiles in this region would have inflections or loops. At s = 0 (Maxwell equation) and 0 < A < 1 there exists a multivalued solution, only.
In order to compute the velocity profile it is necessary to calculate the following integral:
d£
w(q) = | qd£ =J q—dq + C .
T, • • J , d£ 2 (1 -q2 + 3sq2 +sq4) Bearing in mind that from (20) — = ———-^—r-i—1
dq G (1 + q2 )2 analytical expression for the longitudinal velocity in a parametric from
q
(21)
we arrive at the
0
0.00 0.50 1.00 1.50 2.00
G£
Fig. 1. Graphical representation of relationship (20) for different values of s .
4
2
m=
2q(l + sq2)
1 + q2 G
w(q) = -J-s(l + q2 )+(l -s)
ln
1 + q2
+
1 + q2
+c > Gl
(23)
An arbitrary function C is determined from the boundary condition w = 0 and q = qw at £ = 1
C = s(l + ql )-(l-s)
ln
1 + qw2
+
2
1 + qw2
where the dimensionless pressure drop and velocity gradient at the wall are related by the following equation
2
2q 1+sqw = G
2qw-:r = G •
1+qw
The velocity profile can be finally expressed as
q(1 + sq2 )(1 + ql) (1 + q 2 ) (1 + sqw )'
S(q) =
(24)
(25)
qw
(1 + ql)
w(q) =
s (ql - q2)+ (1 -s)
f
ln
1 + q2
1+ql
+
2qw(1 + sql)
2(ql - q2 ) ^ (1 + q2 )(1 + ql)
The obtained velocity profile is then used to define the dimensionless flow rate
1 qw df
Q(q) = 2tcJwf df =2rc J w^ f dq 0 0 dq
and dimensionless stress power
q S AS n f^ff2 At n W ^ ff2 d£ 2
N(q) = Tq£ d£ =n G jd£ =% G dq,
(26)
(27)
(28)
where T = x,
A
Q = 4h£, N = N<^
Mo L
M ro
The dependencies of the flow rate and stress power on the dimensionless pressure gradient are also represented in a parametric form by adding relationship (24) to equations (27) and (28).
1
2
Discussion
In the previous section, we have obtained an analytical solution in a parametric form, which is multivalued in the sense that it describes multitude of transverse velocity profiles including those that lead to tangential discontinuities.
As an example, consider plots of the dimensionless pressure drop G, flow rate Q and stress power N versus parameter qw having the meaning of dimensionless velocity gradient at the capillary wall. Figure 2 shows curves obtained from equations (24), (27) and (28) at s = 0.05. It is seen that there is a range of parameters where the flow rate and stress power fall with increase of qw, which corresponds to unstable unrealistic motion. Moreover, there exists such a range of parameters where the flow rate and stress power are negative. The
1.2 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4
4
Q
Fig. 2. Dimensionless pressure drop G, flow rate Q and stress power N versus dimensionless velocity
gradient at the channel wall qw .
1.2 1.0 w 0.8 0.6 0.4 0.2
0.2
0.4 0.6
£
0.8
1.0
Fig. 3. Velocity profiles at G=Gmax=1.0559 and £ =0.05:
1 - velocity at point A (0=0.4623, #=0.4886);
2 - velocity at point B (0=2.1222, #=2.2409).
velocity profiles corresponding to the discending (A-F) and ascending (F-C) parts of the curve G - qw (dashed line in Fig. 2) are kinematically inadmissible. Figure 4 gives an example of such kinematically inadmissible velocity profile at point E (curve 3). Thus, the realistic motions are those which correspond to the ascending parts of the curve G - qw shown in Fig. 2 by solid lines.
Let us assume that the flow is generated by a prescribed pressure gradient. From Fig. 2 it can be concluded that at certain values of the pressure gradient the problem has two linearly stable solutions and one unstable solution. To answer the question as to which of the stable solutions is realized in nature we must turn to the flow prehistory.
A monotonic growth of the pressure gradient G results in a subcritical bifurcation -the multivaluedness of solutions appears before the pressure gradient reaches the critical value. Since the branched subcritical solutions are unstable at infinitesimal amplitudes of disturbances, the continuos bifurcation is not realized in nature. Instead, one can anticipate the occurrence of discontinuous process, in which the perturbation solution leaving the basin of
w
Fig. 4. Velocity profiles at G=Geq= 0.9111 and £ =0.05:
1 - velocity at point D (0=0.3563, N=0.3247);
2 - velocity at point C (0= 0.3563, N=0.3247);
3 - velocity at point E (0=-0.0007, N=-0.0007).
the basic solution attraction will jump over the unstable solution to a stable one with a higher flow rate and stress power. The velocity profiles for this new linearly stable regimes have either slight tangential discontinuities (continuity of a function and a jump in the function derivative) or drastic tangential discontinuities (discontinuity of a function). The drastic discontinuities are unstable even against infinitely small disturbances [10] so that only slight tangential discontinuities are realizable. The slight discontinuities correspond to a transition along the line CD, which at intersection with the curve G - qw generates two regions of
equal area (Si = S 2).
At monotonic decrease in the pressure gradient an inverse transition to the left branch of the curve G - qw along the line CD occurs at Geq < Gmax. At this pressure drop two stable
solutions yield equal values of the flow rate and stress power. In the second solution, the points of slight discontinuities on the velocity profile are located on the channel boundary (curve 2 in Fig. 4). Inside the channel both solutions are identical. Thus, under condition of monotonic pressure decrease the velocity rearrangement does not occur and the flow rate and stress power do not suffer any jumps.
Since the flow rate jumps due to increase or decrease in the pressure drop take place at different critical values of G the flow rate - pressure drop relationship has a hysteresis (Fig.
5).
G
1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0
0.02 0.04 0.06 0.08 0.10 0.12 0.14
8
Fig. 6. Phase plane G - 8 .
1 - stability region: solution is unique and velocity profile is smooth and nearly parabolic;
2 - metastability region: two linearly stable solutions, which are unstable for finite amplitude disturbances;
3 - stability region: solution is unique and velocity profile has slight tangential discontinuities.
For 8 < 1/9 there is a certain region Geq ^ Gmax, in which the basic motion is
metastable for infinitely small disturbances and unstable for finite amplitude disturbances. The critical condition for flows of this type is defined by the empirical expression
E/^1 - (a + b)(a + c) = 1 + 0.2376 —2, where the third elastic number El3 is expressed
—
as
El3 = — = 1—1——ll. in particular, the critical condition for the DeVitta model is
x e L |
merely El3 = 1. In a subcritical mode, the velocity profiles are continuous, smooth and nearly parabolic, while in a supercritical mode at (a + b)(a + c) < 1 the velocity profiles in the general case have slight tangential discontinuities. For the Maxwell type models (— 2 = 0) at El3^1 - (a + b)(a + c) > 1 kinematically admissible supercritical solutions are nonexistent.
At arbitrary values of s the phase plane G -8 involves three regions of solution: a) the region where the problem has unique stable solution (unshaded region 1 in Fig. 6); b) the region of metastable solution (solution is stable for infinitely small disturbances and unstable for finite amplitude disturbances - shaded region 2 in Fig. 6), where realization of one of the stable solutions is governed by the flow prehistory and c) the region of unique stable solution with slight tangential discontinuities (region 3 in Fig. 6) bounded above by an ordinary hydrodynamic stability limit.
The obtained exact analytical solution predicts some effects, which are characteristic of the blood flow in a vessel. At certain values of the model parameters and at pressure gradient exceeding the critical value the solution predicts a simplified velocity profile, which is common to a near-wall layer with shear rates essentially higher than in the center of the flow. If we suppose that the diffusion motion of erythrocytes is directed from the region of higher shear rates toward that of lower shear rates then we can readily explain the existence of the layer with reduced erythrocyte concentration experimentally observed at the vessel walls. The thickness of this layer can be defined from the obtained solution. The predicted hysteresis effect in the pressure drop-flow rate characteristic is expected to be the reason for doubling the frequency of photometric signal (two-hump effect) in measuring the intensity of light reflected from blood.
1
References
1. ЛЕВТОВ В.А., РЕГИРЕР С.А., ШАДРИНА Н.Х. Агрегация и диффузия эритроцитов. Современные проблемы биомеханики, 9: 5-33, 1994 (in Russian).
2. HINCH E.J., LEAL L.G. Constitutive equations in suspension mechanics. Part 1. General formulation. J Fluid Mech, 71(3): 481-495, 1975.
3. TRUESDELL C., NOLL W. The non-linear fluid theories of mechanics. III/3, Hunbuch der Physics, Berlin-Heidelberg-New York, Springer, 1965.
4. COLEMANN B.D. Kinematic concepts with applications in the mechanics and thermodynamics of incompressible viscoelastic fluids. Arch Rat Mech Anal, 9: 273, 1962.
5. DUNN J.E., FOSDICK R.L. Thermodynamics, stability and boundedness of fluids of complexity 2 and of second grade, Arch Rat Mech Anal, 56: 191, 1974.
6. АНДРЕЙЧЕНКО Ю.А., БРУТЯН М.А., ОБРАЗЦОВ И.Ф., ЯНОВСКИЙ Ю.Г. Спурт-эффект для вязкоупругих жидкостей в 4-константной модели Олдройда. ДАН, 32(3): 327-330, 1997 (in Russian).
7. ULTMANN J.S., DENN M.M. 40th Meeting of the Society of Rheology, St Paul, 1969.
8. HAN CH.D. Reology in polymer processineg. New York-San Francisco-London, Academic Press, 1976.
9. ASTARITA G., MARRUCCI G. Principles of Non-Newtonian fluid mechanics. McGraw-Hill, 1974.
10. ЛАНДАУ Л.Д., ЛИФШИЦ Е.М. Гидродинамика, Москва, Наука, 1986 (in Russian).
11. PETRIE C.J.S., DENN M.M. Instabilities in polymer processing. AIChEJ, 22: 209, 1976.
12. JOSEPH D.D. Stability of fluid motions. Berlin - Heidelberg - New York, Springer - Verlag, 1976.
13. PEARSON J.R.A. Instability in non-Newtonian flow. In: Ann Rev Fluid Mech, Vol. 8, Polo Alto: Annual Reviews Inc, 1976.
ВЯЗКОУПРУГИЕ ЭФФЕКТЫ ПРИ ТЕЧЕНИИ КРОВИ В НЕДЕФОРМИРУЕМОМ КАПИЛЛЯРЕ
С.Н. Аристов, О.И. Скульский (Пермь, Россия)
Полидисперсная структура крови, основу которой составляют эритроциты, приводит к её сложному реологическому поведению. В первом приближении течение крови по сосудам моделируется стационарным движением вязкоупругой жидкости в капилляре с недеформируемыми стенками. Получено точное аналитическое решение пуазейлевского течения вязкоупругой жидкости, поведение которой описывается обобщенной 6-ти константной моделью Джеффри. Профили скорости получены в параметрическом виде, где роль параметра играет градиент скорости. Определены критические значения безразмерного градиента давления, превышение которых приводит к слабым тангенциальным разрывам в продольном профиле скорости. При плавном повышении и понижении градиента давления в некоторой области параметров предсказывается явление гистерезиса. Библ. 13.
Ключевые слова: течение крови, неньютоновская жидкость, точное решение, тангенциальные разрывы, гистерезис
Received 19 November 1999
