The mathematical modelling and the numerical simulations have become important tools for better understanding of the human cardiovascular system in recent years. One of the main goals in investigating the flow in the aortic system is to understand arteriosclerosis and the related phenomena as well as their dependence on a blood flow structure. In particular, the aorta and arteries have a low resistance to blood flow compared with the arterioles and capillaries. When the ventricle contracts, a volume of blood is rapidly ejected into the arterial vessels. Since the outflow to the arteriole is relatively slow because of their high resistance to flow, the arteries are inflated to accommodate the extra blood volume. During diastole, the elastic recoil of the arteries forces the blood out into the arterioles. Thus, the elastic properties of the arteries help to convert the pulsatile flow of blood from the heart into a more continuous flow through the rest of the circulation. Hemodynamics is a term used to describe the mechanisms that affect the dynamics of blood circulation    . In reality, an accurate model of blood flow in the arteries would include the following realistic features: 1) the flow is pulsatile, with a time history containing major frequency components up to the eighth harmonic of the heart period; 2) the arteries are elastic and tapered tubes; 3) the geometry of the arteries is complex and includes tapered, curved, and branching tubes and 4) in small arteries, the viscosity depends upon vessel radius and shear rate  . Such a complex model has never been accomplished. But each of the features above has been “isolated,” and qualitative if not quantitative models have been derived. As is so often the case in the engineering analysis of a complex system, the model derived is a function of the major phenomena one wishes to illustrate.
Our goal is to model and examine the general trend of possible solutions associated with the governing equations describing a simple one-dimensional blood flow that would depict a blood progressing within a thin and elastic pulsating artery. In reality, for many flow situations, the changes of density due to changes in pressure associated with the flow are very small but not zero. In our simulations, the density is assumed variable for the following reason: hereafter we treat blood not as a homogeneous fluid but a suspension of particles (red cells, white cells, platelets) in fluid called plasma. Blood particles must be taken into account in the rheological model in smaller arterioles and capillaries since their size becomes comparable to that of the vessel     . In particular, as has been discussed in  , red blood cells (RBCs) exhibit a unique deformability, which enables them to change shape reversibly in response to an external force. Human RBCs have the ability to undergo large deformations when subjected to external stresses, which allows them to pass through capillaries that are narrower than the diameter of a resting RBC. In fact, RBCs are more deformable than any other biomaterial. RBCs are biconcave discs, typically 6 - 8 μm in diameter and 2 μm thick, and their deformation can involve a change in cell curvature, a uniaxial deformation, or an area expansion. In mammals, RBCs are non-nucleated and consist of a concentrated hemoglobin solution enveloped by a highly flexible membrane. The deformability of RBCs plays an important role in their main function, the transport of gases (O2 and CO2) via blood circulation (see also  and  ).
To put the deformability of blood due to pressure in perspective, consider a multi-component system of total volume V, with
where is the subvolume of component i in the system. The (isothermal) compressibility of the system is
But the compressibility of each component is
Therefore, (2) reduces to
Finally, denoting the volume fraction of each component of the system by , we have
In our case, the main contribution to the compressibility and deformability of the blood is coming from RBCs.
The method of obtaining general solution of the governing equations for the given model is considered from two prospective points of view. The first approach represents a singular perturbation theory providing explicit approximate solutions to the model and the second one is based on group theoretical modeling that provides the exact solutions written in an implicit form. The main goal is to compare these two approaches and fetch out the efficiency and deficiency of each proposed approach.
The range of developed models or models being developed extends from lumped models to complicated three-dimensional fluid-structure models  and  . In this article we consider a simple one-dimensional model of blood flow in a vessel. The blood flow in the vessel is described by this and generally by all one-dimensional models is not suitable for describing blood flow in complicated morphological regions as stenosis or bifurcations. However, these situations can also be covered to certain extent and, from one hand, can be used as an alternative to the more complex three-dimensional fluid-structure models or in conjunction with them in a geometrical multiscale fashion, as explained in  . On the other hand computational complexity of one-dimensional models is several orders of magnitude lower than that of multi-dimensional models. Few decades ago, a multi scale approach has attracted wider interest. Namely, in a multiscale approach, one-dimensional models may be coupled on the one hand with lumped-parameter models  based on a system of ordinary differential equations   , or to three-dimensional fluid-structure models, as discussed in  and  . In the latter case they may also provide a way of implementing more realistic boundary conditions for 3D calculations; or, they can be used for the numerical acceleration of a three-dimensional Navier-Stokes solver in a multilevel-multiscale scheme. Additionally, one-dimensional models give a good description of the propagation of pressure waves in arteries   , hence they can be successfully used to investigate the effects on pulse waves of the geometrical and mechanical arterial modification, due e.g. to the presence of stenoses, or to the placement of stents or prostheses  .
In order to describe a problem in mathematical terms, one must make use of the basic laws that govern the elements of the problem. Within the frame of the present modeling, we start with conservation laws for mass and momentum and consider a perfect compressible fluid propagating along a tube with longitudinal coordinate x and slowly varying cross-section . Because of the pressure gradient in the blood, the artery wall must deform. The elastic restoring force in the wall makes it possible for waves to propagate. In terms of one-dimensional modeling, we assume that the artery radius varies from the constant mean in time and along the artery (in x). We denote the local cross sectional area be , and the averaged velocity be . Consider a fixed geometrical volume between x and , through which fluid moves in and out. Conservation of mass requires
We next assume that the time rate of momentum change in the volume is balanced by the net influx of momentum through the two ends and the pressure force acting on all sides. The rate of momentum change M is given by
where is the density of fluid associated with the mixture density of the blood consisting of blood plasma and red blood cells (RBC). In reality, RBC fraction may include large viscosity variations, stressing the importance of accounting for the non-Newtonian effects (see e.g.  , where the Quemada viscosity model  is used to account for the non-Newtonian viscosity behavior).
The net rate of momentum influx is
The net pressure force at the two ends is given by while that on the sloping wall is
The sum of all pressure forces P is given by
Balancing the momentum by equating M given by (7) to the sum of (8) and P given by (10) we get, after making use of mass conservation (6),
Arterial pulse propagation varies along the circulatory system as a result of the complex geometry and nonuniform structure of the arteries. In order to learn the basic facts of arterial pulse characteristics, we assume an idealized case of an infinitely long circular elastic tube that contains a slightly compressible blood, which is a suspension of particles in what’s basically water. As such, it’s compressibility will be mainly due to the RBC, as explained above. We can think of it as a two-phase homogeneous nonviscous fluid flow of water and gas bubbles. If we apply pressure to the water/gas mixture the overall density will decrease as the gas compresses, leading to the mixture continuity equation that, under the assumption of zero relative velocity, reduces the equivalent single phase flow of density  :
In addition, empirical constitutive laws are needed to relate pressure and density such as equations of state
where S denotes the entropy. Since in our modeling no temperature gradient is imposed externally and the gradient of the flow is not too large, we ignore thermal diffusion. The fluid motion is the adiabatic; entropy is constant. As a result, depends only on density. As we have from thermodynamics,
where T is the temperature and is the ratio of specific heats at constant pressure and constant volume. Furthermore, since pressure is a function of density only, we can write . Expanding this function in a Taylor series about the equilibrium density , we have
where is the equilibrium pressure at which . Since is small, we can neglect the second- and higher-order terms and write
where is a constant. From this equation it follows that
Since the force due to gravity is neglected, combining (11), (12) and (17), we arrive at the governing equations of motion for unknowns velocity , pressure and density that are written as follows:
in which t is time, x is a spatial variable and is a constant.
We eliminate the pressure from these equations by differentiating the Equation (18) with respect to t and using the equation of state (20) to get
Using Equation (19), we can rewrite (21) as
As follows from Equations (18) and (20), we have
Substituting this result into (22), we arrive at the following single equation for :
3. First Approach: Approximate Analysis
In order to identify the resonant input in the model, we start with a an approximate solution in the form of naive expansion
where is small parameter and is an exact trivial solution of Equation (24).
3.1. Failure of the Direct Approach
We substitute the expansion (25) into (24) and collect powers of .
Problem gives the following equation:
We look for solution of Equation (26) in the form
where are constants and in which and are wave number and frequency of the primary wave. Solution (27) is valid provided the dispersion relation is satisfied:
So the problems represent a hyperbolic model with two wave modes
Problem gives the following equation:
In view of presentation (27), the right hand side in Equation (30) is written as
where we denote
We look for particular solution in the form
Substituting this solution into (31), we obtain the following expressions for H and R:
As expected, because of the dispersion relation (28), the right hand side of Equation (30) corresponds to resonance and yields the secular terms.
In particular, if we look for particular solution of the form
the resonance input would be removed since, in this case, the constants and would be
and so the particular solution would have the form
Since grows linearly in time, the term would become comparable
to for large values of time t (e.g. when time is of order ), as shown in
In particular, Figure 1 shows the qualitative behavior of the time series of the second order approximation of solution with secular terms (37) for the values and 2 and and 2. The following values of parameters have been chosen: and is determined by the dispersion relation (29). Since we are interested only in general solutions, the choice of constants is arbitrary and we are focused on qualitative analysis.
As seen from (37), the first terms of the expansion (25) provide a local (small t) approximation, at most. The shortcoming of (25) is related to the breakdown of the straightforward approach on nonlinear perturbation analysis of equation (24), but is more transparent to explanations. The nonlinear terms in (24) will slowly, but accumulative, absorb energy and damp the motion. Hence, even
Figure 1. Visualization of the approximate solution with secular terms in .
though the term itself is small the long term effects are crucial and the solution cannot be described as being periodic plus a small correction. The consequence for a naive expansion (25) is that the ordering requirement is violated. However, it may be instructive to try and fail in order to understand the nature of the resonance phenomena.
3.2. Multiple Scale Approach
We introduce the latter scale according to the new variable
We now consider the fast scale t, and the slow scale , as independent variables. We rewrite Equation (24) in terms of the new variable (38) and modify the series expansion (25) to the form
which yields the perturbation hierarchy similar to (26) and (30), i.e.
where F represents the right hand side if Equation (30) and
appears because of scaling (38). Since the derivative with respect to the fast variable appears only at , the problem at is identical to Equation (26). The slow time variable is implicit in the constant of integration and the most general real-valued solution of the problem can be written as
With this solution in hand, Equation (40) reads
To avoid secular terms we must require
The natural choice is to set Then, squaring and adding the resulting equations for , we arrive at a single equation
Solving Equation (46), we write the solution in the form
where is a constant of integration and is related to and by the dispersion relation (28), i.e.
For the purpose of visualization, Figure 2 is used to compare the qualitative
Figure 2. Visualization of the approximate solution .
behavior of the time series of the of solution (47) for the values of and when and for the same values of parameters and as we used above to visualize the time series of shown in Figure 1.
A question of particular interest is the investigation of asymptotic stability of the trivial solution . This will be done in the next section.
3.3. Stability of Perturbed Steady Flows
We note that the stationary solution
solves Equation (24) in the stationary case, i.e.
which can be integrated to give the exact solution of the form
where and are constants. We denote
Then the only real branch of the solution for of Equation (49) can be written explicitly as
Figure 3 is used to visualize the stationary flow given by (52) for different values of the constant and the following fixed values of parameters: and
Figure 3. Visualization of the stationary solution .
Since, as Figure 2 shows, is growing linearly in x, we classify as non-physical solution.
Let us now look for a nonstationary solution of Equation (24) that is close to in the form (see e.g.   )
where is a small parameter and denotes the perturbation. This procedure is largely formal. Mathematics ideal requires proof that the solution of the complete equations in question for has a solution of the approximate equations at zeroth order of (at least asymptotically). In fact, this ideal is achieved in very rare cases; researchers are usually limited to the formal construction of an approximate model. The justification is based on physical intuition which opens a wide scope. It is clear that, at the same time, the role of the criterion of practice is greatly increased.
We assume a perturbation of the form of a plane harmonic wave propagating in the positive x direction,
where A is a constant amplitude, k is a wave number and is the angular frequency of the oscillations. Substituting the presentations (53) and (54) into Equation (24) and collecting the terms of the order , we get the nonlinear equation for the mean flow
which coincides with (49) and thus has the form (52).
Similarly, collecting and separating the real and imaginary parts of the terms of the order , we get the equations
For progressing wave like solution, Equations (56) and (57) implies that there is another particular exact solution of Equation (49) (and, correspondingly, of Equation (24))
provided that and k satisfy the dispersion relation (28), i.e.
with two known wave modes given by (29). As one can expect, since the flow is away of frictional boundaries, the dispersion relation (28) confirms asymptotic stability of the mean constant flow (58).
Figure 4 shows snap-shots of the perturbed flow at
Figure 4. Visualization of the perturbed constant flow for different values of wave number k.
initial time (red line) and later times of and 1.8 units (plotted in different colors) for different values of the wave number k and the following fixed values of parameters: (for ), (for ), (for ), and (for ), and .
4. Second Approach: Group Theoretical Point of View
Detailed presentations of the theory of symmetries and invariant solutions of differential equations can be found elsewhere     . For convenience, we summarize the basic notation from calculus of Lie group analysis in Appendix, which represents a simplified version of the overview of basic concepts of Lie symmetry groups.
A simple inspection shows that Equation (24) admits the one-parameter groups of translations
of t and x and the one-parameter group of uniform scaling transformations in the -plane:
The above transformations groups have the following generators (called also infinitesimal symmetries):
Various linear combinations of the generators (60) can serve for constructing group-invariant solutions of Equation (60). A brief description of group invariant solutions and illustrative examples are given in  .
4.1. Traveling Waves
Traveling waves are invariant solutions with respect to any linear combination of the translation generators and . We take the linear combination in the form
The operator (61) has two independent invariants, v and According to the theory of invariant solution (see, e.g.  , Section 7.2), the invariant solution has the representation
The representation (62) reduces Equation (24) to the ordinary differential equation
We integrate it once and obtain
The second integration gives the cubic equation
for determining . Thus, the traveling wave solution (62) is determined by the cubic Equation (63) and involves three arbitrary parameters . Cardan’s solution for the cubic equation (see, e.g.  , Section 1.1.1) allows to express the traveling waves in radicals.
Remark: In the special case in (61), Equations (62) and (63) give the stationary solution given explicitly by the cubic Equation (50).
4.2. Similarity Solution
The invariant solution with respect to the generator
of the uniform scaling transformation group is known as a similarity solution. The operator (64) has two independent invariants, v and
The invariant solution has the representation
Calculations show that the representation (66) reduces Equation (24) to the second-order ordinary differential equation
Rewriting this equation in the form
and integrating once, we arrive at the first-order equation
We transform Equation (68) into separable form
by introducing the new dependent variable
Separating the variables,
we write Equation (69) as
whence upon integration
Evaluating the integral in (71) we obtain the following cases.
Case 1: Then the integration in Equation (71) gives
Case 2: Then we evaluate the integral by writing
and write Equation (71) as
Case 3: In this case Equation (71) becomes
Based on the above calculations we conclude that all similarity solutions based on the infinitesimal symmetry (64) are provided by Equations (65), (66), (70) and (72)-(75).
We have investigated a nonlinear partial differential equation of second order that could be used to model a simple one-dimensional blood flow inside a tube of varying cross-section. This model can be an approximation for a pulsating elastic artery. We have proposed two different points of view. The first approach represents a singular perturbation theory that formalizes the scale-separation property by explicitly defining multiple scales that exist in the given nonlinear model with the goal of separating derivatives with respect to fast and slow scales into different orders of perturbation theory. The advantage of this approach is that it yields a solvable perturbative hierarchy of equations that provides useful perturbative information already at low orders in . However, the disadvantage of this approach is the need to identify the various scales a priori and, in the frame of the present modeling, multiple scale approach cannot be brought beyond the leading order. Alternatively, group theoretical approach provides all possible exact solutions of the nonlinear model (24) without any perturbations and, consequently, without introducing a small parameter , which is a significant advantage. In this article, we have provided the exact solutions that were obtained implicitly by solving the nonlinear ordinary differential equations, which have the deficiency of the latter approach. However, in terms of numerical simulations, the second approach seems more advantageous.
This research was supported in part by a Summer Faculty Fellowship Grant provided in Spring, 2017 by College of Natural and Health Sciences, University of Wisconsin-Parkside.
Appendix: Outline of Methods from Lie Group Analysis
Basic concepts from Lie group analysis of differential equations that are used in the present paper are assembled here. For further information regarding Lie groups and their applications to the theory of differential equations, the reader should consult the various classical and modern texts in the field, such as         . A concise introduction to the calculus of Lie symmetry groups can be found in  . Readers interested in applications of Lie groups in fluid dynamics can find a wealth of information in  and  .
Definition of one-parameter groups. Let
be a one-parameter family of invertible transformations of points into points Here a is a real para- meter from a neighborhood of and we impose the condition that Trans- formation (A1) is an identity if and only if i.e.,
The set G of transformations (A1) satisfying Condition (A2) is called a (local) one-parameter group of transformations in if the successive action of two transformations is identical to the action of a third transformation from G, i.e., if the function satisfies the following group property:
with a smooth function defined for sufficiently small a and b. The group parameter a in the transformation (A1) can be changed so that the function (A4) becomes In other words, the group property (A3) can be written, upon choosing an appropriate parameter a (called a canonical para- meter) in the form
Group Generator. Let G be a group of transformations (A1) satisfying the condition (A2) and the group property (A5). Expanding the functions into Taylor series near and keeping only the linear terms in a, one obtains the infinitesimal transformation of the group G:
The first-order linear differential operator
is known as the generator of the group G.
Invariants. A function is said to be an invariant of the group G if for each point is is constant along the trajectory determined by the totality of transformed points
The function is an invariant of the group G with Generator (A8) if and only if
Hence any one-parameter group has exactly functionally independent invariants (basis of invariants). One can take them to be the left-hand sides of first integrals of the characteristic equations for linear partial differential Equation (A9). Then any other invariant is a function of
Invariant equations. We say that a system of equations
is invariant with respect to the group G (or admits the group G) if the trans- formations (A1) of the group G map any solution of Equations (A10) into a solution of the same equations, i.e.,
whenever z solves Equations (A10). The group G with the generator (A8) is admitted by Equations (A10) if and only if
where the symbol means evaluated on the solutions of Equations (A10).
If z is a collection of independent variables dependent variables and partial derivatives of u with respect to x up to certain order, where
then (A10) is a system of partial differential equations
Furthermore, if the transformations (A1) are obtained by the transformations of the independent and dependent variables
and the extension of (A14) to all derivatives etc. involved in the differential Equations (A13), then Equations (A11) define a group G of transformations (A14) admitted by the differential Equations (A13). In other words, an admitted group does not change the form of the system of differential Equations (A13). The generator of the admitted group G is termed an infinitesimal symmetry (or simply symmetry) of the differential Equations (A13). Equations (A12) serve for obtaining the infinitesimal symmetries and are known as the determining equations. These equations are linear and homogeneous and therefore the set L of its solutions is a vector space. Integration of determining equations often provides several linearly independent infinitesimal symmetries. Moreover, the determining equations have a specific property that guarantees that the set L is closed with respect to the commutator Due to this property, L is called a Lie algebra. If the dimension of the vector space L is equal to r, the space is denoted by and is called an r-dimensional Lie algebra. An r-dimensional Lie algebra generates a group depending on r parameters which is called an r-parameter group.
Invariant solutions. Let the differential Equations (A13) admit a multi- parameter group G, and let H be a subgroup of G. A solution
of Equations (A13) is called an H-invariant solution (termed for brevity an invariant solution) if Equations (A15) are invariant with respect to the subgroup H. If H is a one-parameter group and has the generator X, then the H-invariant solutions are constructed by calculating a basis of invariants .
Submit or recommend next manuscript to SCIRP and we will provide best service for you:
Accepting pre-submission inquiries through Email, Facebook, LinkedIn, Twitter, etc.
A wide selection of journals (inclusive of 9 subjects, more than 200 journals)
Providing 24-hour high-quality service
User-friendly online submission system
Fair and swift peer-review system
Efficient typesetting and proofreading procedure
Display of the result of downloads and visits, as well as the number of cited articles
Maximum dissemination of your research work
Submit your manuscript at: http://papersubmission.scirp.org/
Or contact email@example.com