Approximate and Invariant Solutions of a Mathematical Model Describing a Simple One-Dimensional Blood Flow of Variable Density

We examine governing equations that could be used to model a one-dimensional blood flow within a pulsating elastic artery that is represented by a tube of varying cross-section. The model is considered from two perspectives. The first represents a singular perturbation theory providing explicit approximate solutions to the model, and the second is based on group theoretical modeling that provides exact solutions in implicit form. The main goal is to compare these two approaches and lay out the advantages and disadvantages of each approach.


Introduction
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 de-sistance 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 [1] [2] [3].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 [4].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 [5] [6] [7] [8].In particular, as has been discussed in [9], 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 (O 2 and CO 2 ) via blood circulation (see also [10] and [9]).
To put the deformability of blood due to pressure in perspective, consider a multi-component system of total volume V, with where i V 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 i α , 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 [11] and [12].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 [13].
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 [12] based on a system of ordinary differential equations [11] [14], or to three-dimensional fluid-structure models, as discussed in [15] and [16].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 [17] [18], 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 [13].

Modeling
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 ( ) , a x t .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 ( ) , r x t varies from the constant mean 0 r in time and along the artery (in x).We denote the local cross sectional area be 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.[19], where the Quemada viscosity model [20] 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 ( ) 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 ρ [8]: 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 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 0 ρ , we have where 0 p is the equilibrium pressure at which 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 ( ) , v x t , pressure ( ) , p x t and density ( ) that are written as follows: , p λρ = (20) 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 ( )

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 0 cont v = is an exact trivial solution of Equation (24).

Failure of the Direct Approach
We substitute the expansion (25) into (24) and collect powers of ε .
Problem ( ) 1 0 ε gives the following equation: We look for solution of Equation ( 26) in the form ( ) where So the ( ) 0 ε 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 H and R would be ( ) and so the particular solution would have the form 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 1 k = and 2 and 1 λ = and 2. The following values of parameters have been chosen: and 1 ω 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

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.
Problem ( ) 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 ( ) 2 0 ε , the problem at ( ) 26).The slow time variable τ is implicit in the constant of integration and the most general real-valued solution of the ( ) 1 0 ε problem can be written as , , e e .
To avoid secular terms we must require ( ) ( ) The natural choice is to set 1 0. A * = Then, squaring and adding the resulting equations for 1 A , we arrive at a single equation Solving Equation ( 46), we write the solution in the form ( ) ( ) where 1 c is a constant of integration and 1 ω is related to 1 k and λ by the dispersion relation (28), i.e.
( ) For the purpose of visualization, Figure 2 is used to compare the qualitative A question of particular interest is the investigation of asymptotic stability of the trivial solution 0 v .This will be done in the next section.

Stability of Perturbed Steady Flows
We note that the solution ( ) solves Equation ( 24) in the stationary case, i.e.
which can be integrated to give the exact solution of the form ( ) ( ) where M and 2 M are constants.We denote ( ) Then the only real branch of the solution for ( ) 49) can be written explicitly as ( ) Since, as Figure 2 shows, ( ) U x is growing linearly in x, we classify ( ) U x as non-physical solution.
Let us now look for a nonstationary solution of Equation ( 24) that is close to ( ) U x in the form (see e.g.[21] [22]) where ε is a small parameter and v denotes the perturbation.This proce- dure is largely formal.Mathematics ideal requires proof that the solution of the complete equations in question for 0 ε → 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, ( , e , where A is a constant amplitude, k is a wave number and ω is the angular fre- quency of the oscillations.Substituting the presentations (53) and (54) into Equation (24) and collecting the terms of the order 0 ε , we get the nonlinear equa- tion for the mean flow which coincides with (49) and thus

( )
U x 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).

Second Approach: Group Theoretical Point of View
Detailed of the theory of symmetries and invariant solutions of differential equations can be found elsewhere [23] [24] [25] [26].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 The above transformations groups have the following generators (called also infinitesimal symmetries):

Traveling Waves
Traveling waves are invariant solutions with respect to any linear combination of the translation generators 1 X and 2 X .We take the linear combination in the form The operator (61) has two independent invariants, v and .z x mt = − According to the theory of invariant solution (see, e.g.[27], 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 ( ) f µ .Thus, the traveling wave solution (62) is determined by the cubic Equation (63) and involves three arbitrary parameters

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 ( ).
and integrating once, we arrive at the first-order equation ( ) We transform Equation (68) into separable form ( )

Conclusion
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 signifi- cant 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.
averaged velocity be ( ) , v x t .Consider a fixed geometrical volume between x and d x x + , through which fluid moves in and out.

Figure 1 .
Figure 1.Visualization of the approximate solution ( ) , v x t with secular terms in 2 p v .

Figure 2 . 1 A
Figure 2. Visualization of the approximate solution ( ) , v x t .

Figure 3 2 ,Figure 3 .
Figure 3 is used to visualize the stationary flow

Figure 4 Figure 4 .
Figure 4 shows snap-shots of the perturbed flow ( ) ( ) 0 , , v x t v v x t ε = + at 60)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[27].
solution for the cubic equation (see, e.g.[27], Section 1.1.1)allows to express the traveling waves in radicals.Remark: In the special case 0 m = in (61), Equations (62) and (63) give the stationary solution ( ) v U x = given explicitly by the cubic Equation (50).
the representation (66) reduces Equation (24) to the second-order ordinary differential equation