1. Introduction
Gidaspow in 1994 [
1] and Gidaspow and Huang in 2009 [
2] employed and studied a model for blood flow based on the kinetic theory which assumes blood as a mixture of plasma and red blood cells (RBCs). This model allows to explain the Fahraeus–Lindqvist effect, the decrease of RBCs viscosity when the diameter of blood vessels decreases. The cause of this effect is the random RBCs kinetic energy per unit mass, known as granular temperature. When this temperature rises near the wall of the vessel, RBCs oscillate and they migrate from the wall to the center of the vessel. This mechanism is called the shear induced diffusion and it is a driving force connected to a gradient of shear. It is observed with the microscope that the smaller the diameter of the blood vessel, the stronger the shear diffusion. This creates a layer of free cells near the wall and consequently a collapse of the viscosity in the narrow vessels. The role of kinetic theory is to provide a theoretical expression for the blood viscosity. The results of Gidaspow’s paper [
2] are in agreement with the experimental data collected for red blood cell concentrations in narrow vessels. Inspired by this work [
2], many researchers have developed models on blood flow to study physiological and patalogical cases. Huang, in her PhD dissertation, applies Gidaspow’s theory to the transport of LDL and HDL lipoproteins and platelets. In [
3], the pulsing flow in a coronary artery is studied using a typical cardiac waveform. The viscosity values of the red blood cells are numerically determined to appreciate the effect and to identify sites where arteriosclerosis could develop. Fedosov et al. [
4] deal with the changes in the shape, mechanics and physical properties of red blood cells, determined by infectious diseases, environmental and hereditary factors. They evaluate how the biorheology of blood flow varies in subjects suffering from malaria and sickle cell anemia and they compare with healthy cases. Felicetti et al. [
5] numerically treat the detection and the treatment of tumors by recognizing tumor markers on the blood cell surfaces. The information is conveyed on one side by nanomachines that have similar size to white blood cells and they can flow in vessels. On the other side, smart probes are used to detect information extrabody about the localization of tumors. Wu et al. [
6] study benchmark problems including a sudden expansion and various driven slots and crevices through the vessels. Their results turned out to be a predictive model for thrombosis in blood-wetted devices. Lopes et al. [
7] treat the flow of blood in the carotid vessels and compare models in which the vessels are idealized with rigid walls and elastic walls. They focus on the properties of the arteries that become more rigid with the appearance of arteriosclerosis. Furthemore, they speculate that the carotid sinus is the area where arteriosclerosis is most likely to occur, given its low values of Time-Averaged Wall Shear Stress. A 3-phase model, consisting of blood, red blood cells and white blood cells, is formulated by Melka et al. [
8]. In particular, they analyze the aorta and its thoracic branches of a child of 8 years, affected by a congenital heart defect, the coarctation of the aorta. Another 3-phase model, this time consisting of plasma, platelets and red blood cells, is explained by Gidaspow and Chandra [
9]. It is observed numerically that platelets and red blood cells move in opposite directions. The high shear on the walls of the vessel causes an increase in granular pressure and temperature that makes the red blood cells oscillate and push them towards the center. On the contrary, the platelets push towards the walls. From these articles, it is clear how expendable the Gispaspow’s theory is in scientific research.
In this paper, we start from this classical model by Gidaspow and Huang [
2] in order to derive a hyperbolic system in the context of Rational Extended Thermodynamics [
10,
11,
12]. Following this non-equilibrium theory, we consider as field variables not only the classical ones in [
2], which are mass densities, velocities and RBCs temperature, but also the stress tensors, the RBG dynamic pressure and the heat flux. In order to describe this extended set of field variables, we need to introduce a new system of field equations. Following the idea of RET, we introduce for these fields a new set of balance laws, and, in particular, here we are inspired by a two-hierarchy structure model recently introduced in [
11,
12] for polyatomic gases. The balance equations needs closure and, following again the main idea of RET, we close them assuming local and instantaneous constitutive relations and using well-known universal physical principles, such as the entropy principle (existence of an entropy inequality and concavity of the entropy density) and the principle of relativity. The exploitation of the principle of relativity helps us in the determination of the velocity-dependence of the constitutive relations, while the principle of entropy implies some useful mathematical relations. The mathematical exploitation of the entropy principle is made by use of the Lagrange multipliers [
13] (main field). The constitutive functions are all determined except for the thermal and caloric equation of state and the production terms.
We identify the remaining terms using the expressions of the RBCs pressure, viscosity and heat conductivity determined by use of the kinetic theory in [
1]. In this way, we assure that our equations coincide with the classical GH model when we are in presence of slowly changing fields and small gradients. In addition, our model can describe other processes where rapid time changes or when a strong deviation from equilibrium occur. In fact, it has been shown that the field equations of RET can describe a range various non-equilibrium phenomena such as light shattering, sound waves, heat waves, structure of shock waves [
10,
11,
12]. RET has been applied to monoatomic gases [
10] and mixtures [
10,
11,
14,
15] with many interesting results. Recently, RET has been generalized to dense and rarefied polyatomic gases both in the classical [
11,
12,
16,
17,
18] and in the relativistic framework [
12,
19,
20], for metal electrons [
21,
22], to quantum systems [
23] and also for biological systems [
24,
25,
26,
27], providing relevant results and good agreement with experimental data.
The model herein introduced, as all ET systems [
28,
29], is symmetric hyperbolic. Hyperbolicity guarantees finite speeds of propagation, while symmetric hyperbolic systems imply well-posedness of Cauchy problems (i.e., existence, uniqueness and continuous dependence on the data). Hyperbolic models overcome the paradox of infinite propagation speed of disturbances, and they are better suited to describe transient regimes. Moreover, they offer a richer scenario of dynamics which can be useful to describe biological problems, see for example [
24,
25,
26,
27] and the references therein. We expect that the model herein introduced can also describe phenomena in blood flow in which a rapid change in the blood pressure occurs, such as the sudden occlusion of a vessel or rapid initial changes in the fields due to external conditions.
Furthermore, from the results obtained for the heat transfer problem in a monoatomic gas or gas-mixtures, see for example [
14,
15,
30,
31] and the references therein, we also expect that our model can show more interesting spatial effects also for time-independent phenomena. Indeed, it has been shown that the field equations of RET are able to put in evidence more spatial effects than classical Navier–Stokes–Fourier equations, when a radial symmetry is present or a particular geometry is taken into account. Therefore, we expect that our model can describe more accurately the blood flow in a vessel and in particular in the case of change in the radial symmetry, like stenosis, or other defects of vessels.
Finally, it was recently shown that the field equations of RET are to describe the thermal-diffusion effect, contrary to the classical equations. In consideration of this RET result, we expect that our model is able to show some interesting additional features to the description of the Fahraeus–Lindqvist effect, which is anyway already present in the Gidaspow classical model [
2].
The remainder of this article is the following: In
Section 2, we recall the Gidaspow classical model, while in
Section 3, we introduce the balance equations for our model.
Section 4,
Section 5 and
Section 6 are devoted to the applications of the physical principles for the determinations of the constitutive relations. In particular, the Galilean invariance principle helps us to determine the dependence of the constitutive relations on the velocity, while the entropy principle implies some algebraic relations for the constitutive relations. In order to avoid too long calculations, we limit the study near the thermodynamic equilibrium by expanding the constitutive relations in terms of the non-equilibrium fluxes. In
Section 7, a summary of the consequences of the principle is made, and in
Section 8, the comparison with the Gidaspow model is shown with the identification of the last unknown functions. So all coefficients and functions here introduced can be obtained in terms of those calculated by Gidaspow using the kinetic theory. In
Section 9, some simple calculations are presented in order to show with a simple numerical example the applicability of our model. Finally, in
Section 10, some conclusions and final remarks are discussed.
2. Gidaspow Model
Gidaspow’s model [
2] describes blood as a mixture of a fluid and the RBCs and uses as field equations the two conservation laws of mass for the two constituents, the two balance laws of momentum and a balance law for the total energy of RBCs. These equations, typical of classical thermodynamics, can be recast in the following form:
where
and
are respectively the plasma and the RBCs densities,
and
their velocities,
and
the stress tensors,
is the RBCs specific interval energy and
is the RBCs heat flux. The field Equation (
1) are closed by the classical Navier–Stokes and Fourier equations, which assume the stress tensor and the heat flux proportional to the gradient of the velocity or the temperature:
In the previous equations and through the paper, the round brackets in the indexes denote the symmetric parts of a tensor. The viscosities
,
and
and the heat conductivity
are calculated by Gidaspow by use of the kinetic theory and written explicitly in [
2]. The right-hand side terms in Equation (
1) are the body forces
due to external specific force
, the actions that the fluid exerts on the RBCs
and vice versa (
is the volume fraction of RBCs). Furthermore, the production terms take into account the interaction between plasma and RBCs, that are
Moreover, the coefficients
and
are evaluated in [
2].
Equations (
1)–(
3) form the classical system of field equations, which consists of nine equations in the nine fields
,
and
.
3. Balance Equations
Starting from Gidaspow’s model, we introduce a system of hyperbolic balance equations. First of all, we consider as field variables not only the densities, the velocities and the temperature but also the stress tensors and the heat flux. The fields equations for all these variables are the two conservation laws of mass for the two constituents, the two balance laws of momentum and the balance law for the total energy of RBCs, as in (
1). Furthermore, instead of assuming the classical Navier–Stokes and Fourier expressions (
2) for the stress tensors and for the heat flux, we introduce for these new variables supplementary balance equations. In particular, the system of balance equations that we assume has the following form:
with
or
s.
Equation (
5)
and (
5)
with
coincide exactly with the classical Gidaspow Equation (
1) indicating with
total energy and
the total energy flux. Equation (
5)
with
is a supplementary equation here introduced for the total energy for the fluid. The other equations of (
5) are introduced assuming that the fluxes in the first equations are the densities in the new equations. This structure was obtained by integration from the Boltzmann equation for a monatomic gasses [
10] and recently this structure was modified for polyatomic gasses [
11,
12,
17], where the two hierarchies structure (for the
F’s and the
M’s) is introduced by considerations from the kinetic theory.
System (
5) represents a set of 24 equations for the 24 field variables that are the two densities
, the two velocities
, the two temperatures
, the traceless (here and in the following, square brackets in the indexes denote the traceless part of a symmetric tensor) part stress tensors
, the RBCs dynamic pressure
and heat flux
. For sake of simplicity, we assume that
. The aim of this paper is to close system (
5), that is to express all quantities in (
5) in terms of these 24 fields.
First of all, in the spirit of Rational Extended Thermodynamics, we assume that the unknown quantities in (
5) are local and instantaneous. This means that is at one time and one point they depend on the field variables at the same time and at the same point. So, indicating a generic unknown quantity with
, it must hold
Then, if the functions
are known, we can insert them into the balance Equation (
5) in order to obtain a closed system of 24 field equations in the 24 fields. The obtained system, for the assumption in (
6), will be a set of quasi-linear first order differential equations, whose solution is called thermodynamic process.
4. Decomposition in Convective and Non Convective Terms
We start now to decompose the densities, fluxes and productions in (
5) into convective and non-convective parts, according to the single component case. Indeed, assuming that the balance equations for each constituent of the mixture must be invariant under a Galilean transformation, it is possible to determine the velocity dependence of the quantities in (
5) on the two velocities, so we have:
The explicit calculations for expressions (
7) are shown in [
10] pp. 53, 54.
It is easy to recognize the decomposition (
7)
of the moment fluxes in the stress tensors and velocity dependent parts, the decomposition (
7)
of the total energies in internal and kinetic energies and the decomposition (
7)
for the total energy flux
in the sum of the heat flux and the velocity dependent parts. These three decompositions are already present in (
1). Then, we may say that
in (
7) are the stress tensors which, according to Gidaspow’s model, are expressed by
while
is the RBCs heat flux.
Then, by insertion of the decomposition (
7) into the balance Equation (
5), one obtains
where the two material derivatives
have been introduced. It can be clearly seen that these equations do not depend on the velocities but only on their derivatives and this is in agreement with the Galileian invariance applied to each constituent of the mixture.
With these decompositions, the unknown functions become
,
,
and
. These quantities must be determined in terms of the set
,
,
,
,
and, in virtue of the principle of material frame indifference applied to the whole mixture ([
10] p. 86, [
11] p. 313, [
12,
32] pp. 6, 177, [
33]), they cannot depend on the velocity of the two components
but only on their relative velocity
. Equivalently, they can be considered as functions of the diffusion velocity
is the mixture velocity. Summarizing, the constitutive functions to be determined are
5. Entropy Principle
The entropy principle asserts that it exists a quantity
h, called entropy, which satisfies a balance equation of the form
for all thermodynamic processes, that is for all solutions of the field equations. In this paper, in line with the methods of Rational Extended Thermodynamics [
10,
11,
12], we use this physical principle in order to determine the constitutive Equation (
11).
First of all, as described before for the fields, the concave entropy density
h, the entropy flux,
, and the entropy production Φ must depend on the field variables through constitutive functions of the form
Then, since the entropy principle must hold for all thermodynamic processes that is for all solutions of the field Equation (
9), we use the field Equation (
9) as constraints for the validity of the entropy principle, This can be mathematically done, by introduction [
13] of the so-called Liu–Lagrange multipliers
’s, so the entropy inequality (
12) becomes
which must be valid for all values of
,
,
,
,
and
.
We start now the mathematical evaluations of (
14). We have to say that the introduction of the Lagrange multipliers adds new unknown quantities. In agreements with the other quantities, we assume that also the
’s must depend on the fields through constitutive relations of the form (
11) and (
13).
It is easy to see that in (
14) there is only one term which depends on the external force
so, since the constitutive relations must be valid independently of the forces applied to the system, we can conclude that
must hold.
Then, in equality (
14) there are only few terms that depend only on the fields and not on their derivatives. These terms are the production terms in the balance equations and the entropy production, so we recover the so-called residual inequality, that yields the entropy production in terms of the other quantities
and whose sign must be verified when the other quantities are determined.
Another relation, that can be easily determined form (
14), is the following:
This relation comes out eliminating all coefficients of the velocity
from equality (
14).
We can now insert the first relations (
15)–(
17) into the balance law (
14) in order to obtain a simpler equation. A very accurate evaluation of the complete entropy principle requires a lot of complex calculations that we prefer not to present here. In the next section, we start by limiting our analysis to processes near the thermodynamic equilibrium. A more general case will be presented in [
34].
6. Evaluation of the Quantities Near Equilibrium
As explained before, from now on we limit our attention to processes close to thermodynamic equilibrium, characterized by vanishing non-equilibrium fluxes
,
,
and
. This can be done by mathematically expanding the constitutive relations around the thermodynamic equilibrium and considering the terms which are at least quadratic in the non-equilibrium fluxes. In this way the entropy quantities become
the constitutive relations are given by
and the Lagrange multipliers take the form
So the determination of the constitutive functions reduces to the calculations of all coefficients in (
18)–(
20), which depends only the equilibrium variables
.
We then proceed by insertion of relations (
15)–(
20) into the entropy inequality (
14) which becomes a very long but linear equation in the derivatives of the field variables
,
,
,
,
and
. As already said, Equation (
14) must be valid for all field variables; therefore, it must be valid also for all gradients of these variables. Then, setting the coefficients of the derivatives of
,
,
,
,
and
equal to zero, we obtain a big set of algebraic equations for the determination of the functions in (
18)–(
20). We prefer not to present all these equations here, but an interested reader can find similar calculations in [
10,
11,
12,
34]. We present here only some equilibrium relations and the results.
6.1. Gibbs Equation
We evaluate the coefficients in expansions (
18)–(
20), starting with the equilibrium terms. Indeed, from the derivatives of
and
in (
18)–(
20) we obtain
which can be recasted in the differential form
This last equation can be compared with the well-known Gibbs relation, that expresses the variation of the equilibrium entropy density in terms of the thermodynamic variables. The expression of the Gibbs relation appropriate to each component of a mixture can be found for example in [
10] p. 90 and it assumes the form
Then, by comparison of (
22) and (
23), taking also into account that
, it is possible to determine the equilibrium part of the Lagrange multipliers, i.e.,
where
represents the chemical potential.
6.2. Evaluation of the Functions Outside Equilibrium
By use of the Lagrange multipliers at equilibrium (
24) it is possible to evaluate from the entropy inequality all coefficients in the constitutive relations (
14) and (
18)–(
20) as explained before. We have the entropy quantities
the first order terms in the constitutive relations
and in the Lagrange multipliers
with
The symbol
represents the specific heat at constant volume, defined as
. We underline that these constitutive relations are in agreement with relations (5.58) p. 123 of Ruggeri and Sugiyama [
11] since the structure for RBCs equations is equivalent to their two hierarchies system. We conclude that all functions introduced in the present analysis are determined except for the pressure
, the internal energy
and the production terms.
6.3. Production Terms
In this subsection, we will limit the generality of the production terms indeed we assume also for the production terms linear constitutive functions and in particular we apply the BGK assumption [
35]. Therefore, the production term in each equation is proportional to the density of the same equation minus its equilibrium value. Explicitly we have:
where
,
and
are the well-known relaxation times that are functions of
and
. In the next section we identify these relaxation times by comparison with the Gidaspow’s model [
2].
8. Comparison with Gidaspow’s Model
In this section, we finally determine the remaining functions by comparison with the Gidaspow–Huang model [
2]. First of all, we identify the RBCs pressure, that by calculations [
1] in kinetic theory is given by
where
is the RBCs volume fraction,
their radial distribution function, that can be chosen according to the Bagnold’s equations, and
e is the restitution coefficient, a measure of the elasticity of the particle collisions (see [
1,
2,
36] for more details). We inserted in
Table 1 some values for these quantities taken in [
2]. Furthermore, we assume that the RBCs energy density is given by
, so that
.
By use of the expression for the RBCs pressure (
30), it is possible to evaluate all functions in (
28) in terms of the RBCs density and temperature, which are
In order to evaluate the relaxation times of the model, we substitute the constitutive relations (
26) and the production terms (
29) into the balance Equation (
9)
and we proceed with the first Maxwellian iterations [
10]. This method is really simple to be used, but it is very useful: In practice, we neglect from the left-hand side of Equation (
9)
the non-equilibrium terms, that is, we consider the case in which no-rapid changes in the fields occur. The obtained equations become
which have exactly the form of the stress tensor, the dynamic pressure and the heat flux predicted by the Navier–Stokes and the Fourier laws (
2). This method is particularly useful; indeed, we can firstly evaluate, by comparison of (
2) with (
32), our relaxation times in terms of the Gidaspow viscosities and heat conductivity. Furthermore, we make sure that our equations predict the same results of the classical model when no-rapid changes in the fields occur. From the comparison and taking into account the expression of the viscosities and heat conductivity obtained in [
2], we easily obtain
where
represents the RBCs diameter, while
represents the fluid volume fraction with
.
Finally, the coefficients
and
in the expression of the interaction terms
and
in (
9), are taken from the analogous terms in [
2] so we set
The first relation takes into account only the linear term of [
2] that is more appropriate for the linear constitutive equations presented in this paper.
All coefficients of our model are evaluated in terms of the field variables; furthermore, with the obtained functions,
is a convex function and the residual inequality (
16) is verified.
Summarizing, system (
9) with the relations (
8), (
26), (
30) and (
31), the productions terms (
29) and (
33), the terms due to the external forces in (
3) and the interaction terms (
4) and (
34) is our system of field equations. It is a closed system of 24 field equations for the 24 field variables
,
,
,
,
and
.
9. A Simple Numerical Example
In order to show clearly the applicability of our model, we present here a simple numerical example. We stress that these results are only a first simple application of our equations and they are not exhaustive of all the applicability of our model. They are also part of another paper [
37], so an interested reader can found there more details.
We use the field equations herein obtained in order to describe the stationary blood flow in a vessel, which is idealized as the gap between two infinite parallel plates. We can assume that the flow is stationary and the fields depend only on the x-coordinate orthogonal to the two walls, so the complete set of field equations becomes a more simple set of ODE. We also assume that the RBCs and fluid velocities have the direction of the flow and they are driven by an external pressure, considered as external force.
It order to recover easily the solution, we linearize the ODE around an equilibrium state characterized by constant
,
and vanishing fluxes. We illustrate in
Figure 1 and
Figure 2 the solution obtained for
, with the parameters in
Table 1, and in terms of the dimentionless variables
The solutions show that the RBCs temperature decreases in the center of the vessel, where instead the RBCs density increases. This result is in agreement with the Fahraeus–Linqvist effect, as already shown in [
2]. The two velocities are illustrated in
Figure 2, showing that the fluid anticipates RBCs, as they are carried from the fluid self.
The solutions in
Figure 1 and
Figure 2 are somehow less general than those presented in [
2], since here only the 1D plane case is considered, no-radial or 2D geometry, nor slip at the boundary as in [
2] and furthermore the equations are linearized around the constant state. Anyway, integration of this very simple and linearized system yields (see for details [
37]) already the presence of some non-vanishing stress components. These non-vanishing components do not affect visibly the classical fields in
Figure 1 and
Figure 2, due to the strong assumptions in this chapter, but they will play surely an important role in all fields when all simplifying assumptions will be removed. This will be shown clearly in further research.
10. Conclusions and Final Remarks
In this paper, starting from a model for blood flow obtained by Gidaspow and Huang [
2] in Classical Thermodynamics, we derived an hyperbolic model based on the methods proper of Rational Extended Thermodynamics [
10,
11,
12]. The stress tensors and the heat flux are considered as new field variables and a new set of field equations is assumed. The introduction of new equations involves new quantities that are here determined as constitutive functions in terms of the fields. This was done using the Galileian and entropy principles. We determine the remaining functions with Gidaspow’s coefficients obtained in [
1]. All parameters herein introduced are determined explicitly in terms of the fields and a simple numerical solution is shown in order to illustrate the applicability of our model.
The model herein obtained is hyperbolic, so it overcomes the paradox of infinite propagation speed of disturbances typical of parabolic ones. Our model can describe transient regimes better than the classical one and it offers a richer scenario of solutions as many other hyperbolic systems. Furthermore, by construction our model converges to the classical one [
2] when the field gradients are negligible. On the contrary, we expect a better description of the phenomena in blood flow when strong changes in the blood occur. As seen for gases and mixture flow, we also expect that our Extended Thermodynamics model can describe more properly situations in with the domain presents more curvature of changes in the radial symmetry.
It must be said that, unfortunately, the integration of our field equations is more complicate than the classical ones, since our equations contain in practice the differential Equation (
9)
instead of the classical Navier–Stokes and Fourier laws (
2). This will imply a more complicate (normally numeric) integration of the whole set of equations and the introduction of initial values for stress tensors and heat flux.
Finally, our model has the same range of applicability of the classical model [
2], so it can describe a great range of blood phenomena and cases. Some applications of our model are already under investigation and they will be part of further works.