JAMP  Vol.5 No.6 , June 2017
Approximate and Invariant Solutions of a Mathematical Model Describing a Simple One-Dimensional Blood Flow of Variable Density
ABSTRACT
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.

1. 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 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 [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 (O2 and CO2) 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

V = i V i (1)

where V i is the subvolume of component i in the system. The (isothermal) compressibility of the system is

κ = 1 V ( V P ) T = 1 V i ( V i P ) T (2)

But the compressibility of each component is

κ i = 1 V i ( V i P ) T (3)

Therefore, (2) reduces to

κ = 1 V i V i κ i (4)

Finally, denoting the volume fraction of each component of the system by α i , we have

κ = i α i κ i (5)

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] .

2. 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 r 0 in time and along the artery (in x). We denote the local cross sectional area be a = π r 2 , and the averaged velocity be v ( x , t ) . Consider a fixed geometrical volume between x and x + d x , through which fluid moves in and out. Conservation of mass requires

a t + ( v a ) t = 0. (6)

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

M = ( ρ v a ) t , (7)

where ρ ( x , t ) 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

( ρ v 2 a ) x d x = ρ v ( v a ) x ρ v a v x . (8)

The net pressure force at the two ends is given by ( p a ) / x while that on the sloping wall is

2 π r p r x = p a x . (9)

The sum of all pressure forces P is given by

P = a p x . (10)

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),

v t + v v x = 1 ρ p x . (11)

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] :

ρ t + x ( ρ v ) = 0. (12)

In addition, empirical constitutive laws are needed to relate pressure and density such as equations of state

p = p ( ρ , S ) , (13)

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 S = S 0 is constant. As a result, p = p ( ρ , S 0 ) depends only on density. As we have from thermodynamics,

( p ρ ) S = γ ( p ρ ) T , (14)

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 p = p ( ρ ) . Expanding this function in a Taylor series about the equilibrium density ρ 0 , we have

p = p 0 + p ( ρ 0 ) 1 ! ( ρ ρ 0 ) + p ( ρ 0 ) 2 ! ( ρ ρ 0 ) 2 + (15)

where p 0 is the equilibrium pressure at which ρ = ρ 0 . Since ρ ρ 0 is small, we can neglect the second- and higher-order terms and write

p = p 0 + λ ( ρ ρ 0 ) (16)

where λ is a constant. From this equation it follows that

p x = λ ρ x and p t = λ ρ t (17)

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 ρ ( x , t ) that are written as follows:

v t + v v x = 1 ρ p x , (18)

ρ t + x ( ρ v ) = 0 , (19)

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

2 v t 2 + v t v x + v 2 v x t = λ ρ ( 2 ρ x t 1 ρ ρ t ρ x ) . (21)

Using Equation (19), we can rewrite (21) as

2 v t 2 + x ( v v t ) = λ 2 v x 2 + x ( λ ρ v ρ x ) . (22)

As follows from Equations (18) and (20), we have

λ v ρ ρ x = v v t v 2 v x . (23)

Substituting this result into (22), we arrive at the following single equation for v ( x , t ) :

2 v t 2 + x ( 2 v v t + v 2 v x ) = λ 2 v x 2 . (24)

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

v ( x , t ) = v 0 + i = 1 ε i v i ( x , t ) , (25)

where ε is small parameter and v 0 = cont 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 0 ( ε 1 ) gives the following equation:

2 v 1 t 2 λ 2 v 1 x 2 + x ( 2 v 0 v 1 t + v 0 2 v 1 x ) = 0. (26)

We look for solution of Equation (26) in the form

v 1 ( x , t ) = A 1 sin θ 1 + B 1 cos θ 1 , (27)

where A 1 , B 1 are constants and θ 1 = k 1 x ω 1 t , in which k 1 and ω 1 are wave number and frequency of the primary wave. Solution (27) is valid provided the dispersion relation is satisfied:

ω 1 2 2 v 0 k ω 1 k 1 2 ( λ v 0 2 ) = 0. (28)

So the 0 ( ε 1 ) problems represent a hyperbolic model with two wave modes

ω 1 = k 1 ( v 0 ± λ ) . (29)

Problem 0 ( ε 2 ) gives the following equation:

2 v 2 t 2 λ 2 v 2 x 2 + x ( 2 v 0 v 1 t + v 0 2 v 2 x ) = x ( 2 v 1 v 1 t + 2 v 0 v 1 v 1 x ) . (30)

In view of presentation (27), the right hand side in Equation (30) is written as

x ( 2 v 1 v 1 t + 2 v 0 v 1 v 1 x ) = δ ( A 1 2 B 1 2 ) cos ( 2 θ 1 ) 2 δ A 1 B 1 sin ( 2 θ 1 ) , (31)

where we denote

δ = 2 k 1 ( ω 1 v 0 k 1 ) . (32)

We look for particular solution in the form

v 2 ( p ) = H cos ( 2 θ 1 ) + R sin ( 2 θ 1 ) . (33)

Substituting this solution into (31), we obtain the following expressions for H and R:

( H , R ) = δ ( B 1 2 A 1 2 , 2 A 1 B 1 ) 4 [ ω 1 2 2 v 0 k ω 1 k 1 2 ( λ v 0 2 ) ] . (34)

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

v 2 ( p ) = t H ˜ cos ( 2 θ 1 ) + t R ˜ sin ( 2 θ 1 ) , (35)

the resonance input would be removed since, in this case, the constants H ˜ and R ˜ would be

( H ˜ , R ˜ ) = k 1 ( A 1 B 1 , A 1 2 B 1 2 2 ) (36)

and so the particular solution would have the form

v 2 ( p ) = k t A 1 B 1 cos ( 2 θ 1 ) k 1 2 t ( A 1 2 B 1 2 ) sin ( 2 θ 1 ) . (37)

Since v 2 ( p ) grows linearly in time, the term ε 2 v 2 would become comparable

to ε v 1 for large values of time t (e.g. when time is of order 1 ε ), as shown in

Figure 1.

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 k = 1 and 2 and λ = 1 and 2. The following values of parameters have been chosen: ε = 0.1 , v 0 = 1 , x = 2 , A 1 = 0.18 , B 1 = 0.19 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

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

though the term ε v 2 ( p ) 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 v 0 > ε v 1 > 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

τ = ε t . (38)

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

v ( x , t , τ ) = v 0 + i = 1 ε i v i ( x , t , τ ) , (39)

which yields the perturbation hierarchy similar to (26) and (30), i.e.

Problem 0 ( ε 1 ) :

2 v 1 t 2 λ 2 v 1 x 2 + x ( 2 v 0 v 1 t + v 0 2 v 1 x ) = 0

and

Problem 0 ( ε 2 ) :

2 v 2 t 2 λ 2 v 2 x 2 + x ( 2 v 0 v 1 t + v 0 2 v 2 x ) = F + Q , (40)

where F represents the right hand side if Equation (30) and

Q = 2 [ 2 v 1 t τ v 0 2 v 1 x t ] (41)

appears because of scaling (38). Since the derivative with respect to the fast variable appears only at 0 ( ε 2 ) , the problem at 0 ( ε 1 ) 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 0 ( ε 1 ) problem can be written as

v 1 ( x , t , τ ) = A 1 ( τ ) e i θ 1 + A 1 e θ 1 . (42)

With this solution in hand, Equation (40) reads

F + Q = 2 ( ω 1 v 0 k 1 ) [ i ( d A 1 d τ e i θ 1 d A 1 d τ e i θ 1 ) 2 k 1 ( A 1 2 e 2 i θ 1 + A 1 2 e 2 i θ 1 ) ] . (43)

To avoid secular terms we must require

( d A 1 d τ + d A 1 d τ ) sin θ 1 = 2 k 1 ( A 1 2 + A 1 2 ) cos ( 2 θ 1 ) , (44)

( d A 1 d τ d A 1 d τ ) cos θ 1 = 2 k 1 ( A 1 2 A 1 2 ) sin ( 2 θ 1 ) . (45)

The natural choice is to set A 1 = 0. Then, squaring and adding the resulting equations for A 1 , we arrive at a single equation

d A 1 d τ = 2 k 1 A 1 2 . (46)

Solving Equation (46), we write the solution in the form

v ( x , t ) = v 0 + ε cos k 1 x ω 1 t 2 k 1 ε t + c 1 + 0 ( ε 2 ) , (47)

where c 1 is a constant of integration and ω 1 is related to k 1 and λ by the dispersion relation (28), i.e.

ω 1 = k 1 ( v 0 ± λ ) .

For the purpose of visualization, Figure 2 is used to compare the qualitative

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

behavior of the time series of the of solution (47) for the values of k 1 = 1 and k 2 = 4 when c 1 = 1.5 and for the same values of parameters ε , v 0 , x and A 1 as we used above to visualize the time series of v 2 ( p ) shown in Figure 1.

A question of particular interest is the investigation of asymptotic stability of the trivial solution v 0 . This will be done in the next section.

3.3. Stability of Perturbed Steady Flows

We note that the stationary solution

v = U ( x ) (48)

solves Equation (24) in the stationary case, i.e.

x ( v 2 v x ) = λ 2 v x 2 . (49)

which can be integrated to give the exact solution of the form

U 3 ( x ) 3 λ U ( x ) = M 1 x + M 2 , (50)

where M 1 and M 2 are constants. We denote

σ = 4 M 1 x + 4 M 2 + 4 ( M 1 x ) 2 + 2 M 1 M 2 x 4 λ 3 + M 2 2 (51)

Then the only real branch of the solution for U ( x ) of Equation (49) can be written explicitly as

U ( x ) = 1 2 σ 2 3 + 4 λ σ 1 3 . (52)

Figure 3 is used to visualize the stationary flow U ( x ) given by (52) for different values of the constant λ = λ 1 , λ = λ 2 and the following fixed values of parameters: λ 1 = 5 , λ 2 = 1 , M 1 = 5 and M 2 = 1.

Figure 3. Visualization of the stationary solution U ( x ) .

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] )

v ( x , t ) = U ( x ) + ε v ˜ ( x , t ) , (53)

where ε is a small parameter and v ˜ denotes the perturbation. This procedure 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,

v ˜ ( x , t ) = A e i ( k x ω t ) , (54)

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 ε 0 , we get the nonlinear equation for the mean flow

x ( U 2 U x ) λ 2 U x 2 = 0 , (55)

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

ω 2 2 U ω k + k 2 ( U λ ) 2 x ( U U x ) = 0 (56)

and

ω k x [ U ( 1 U ) ] = 0. (57)

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))

v ( x , t ) = v 0 = const , (58)

provided that ω and k satisfy the dispersion relation (28), i.e.

ω 2 2 v 0 ω k + ( v 0 λ ) k 2 (59)

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 v ( x , t ) = v 0 + ε v ˜ ( x , t ) at

Figure 4. Visualization of the perturbed constant flow v 0 for different values of wave number k.

initial time t = 0 (red line) and later times of t = 0.5 , 1.0 , 1.5 and 1.8 units (plotted in different colors) for different values of the wave number k and the following fixed values of parameters: v 0 = 1 (for k = 1 ), v 0 = 2 (for k = 2 ), v 0 = 3 (for k = 10 ), and v 0 = 4 (for k = 25 ), a = 1 and ε = 0.1 .

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 [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

t ¯ = t + a 1 , x ¯ = x + a 2

of t and x and the one-parameter group of uniform scaling transformations in the ( t , x ) -plane:

t ¯ = t e a 3 , x ¯ = x e a 3 .

The above transformations groups have the following generators (called also infinitesimal symmetries):

X 1 = t , X 2 = x , X 3 = t t + x x (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] .

4.1. Traveling Waves

Traveling waves are invariant solutions with respect to any linear combination of the translation generators X 1 and X 2 . We take the linear combination in the form

X 1 + m X 2 = t + m x , m = const . (61)

The operator (61) has two independent invariants, v and z = x m t . According to the theory of invariant solution (see, e.g. [27] , Section 7.2), the invariant solution has the representation

v = f ( z ) . (62)

The representation (62) reduces Equation (24) to the ordinary differential equation

( m 2 λ ) f 2 f ( 2 m f f ) = 0.

We integrate it once and obtain

( m 2 λ ) f 2 f 2 m f f = M 1 , M 1 = const .

The second integration gives the cubic equation

f 3 3 m f 2 + 3 ( m 2 λ ) f = M 1 μ + M 2 (63)

for determining f ( μ ) . Thus, the traveling wave solution (62) is determined by the cubic Equation (63) and involves three arbitrary parameters m , M 1 , M 2 . Cardan’s 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 m = 0 in (61), Equations (62) and (63) give the stationary solution v = U ( x ) given explicitly by the cubic Equation (50).

4.2. Similarity Solution

The invariant solution with respect to the generator

X 3 = t t + x x (64)

of the uniform scaling transformation group is known as a similarity solution. The operator (64) has two independent invariants, v and

ξ = x t (65)

The invariant solution has the representation

v = g ( ξ ) . (66)

Calculations show that the representation (66) reduces Equation (24) to the second-order ordinary differential equation

[ ( g μ ) 2 λ ] g + 2 ( g μ ) g 2 2 ( g μ ) g = 0.

Rewriting this equation in the form

[ ( g μ ) 2 g ] λ g = 0 (67)

and integrating once, we arrive at the first-order equation

[ ( g μ ) 2 λ ] g = C 1 , C 1 = const . (68)

We transform Equation (68) into separable form

( w 2 λ ) d w d μ = w 2 + λ + C 1 (69)

by introducing the new dependent variable

w = g μ . (70)

Separating the variables,

w 2 λ λ + C 1 w 2 d w = d μ ,

we write Equation (69) as

d w + C 1 λ + C 1 w 2 d w = d μ ,

whence upon integration

C 1 d w λ + C 1 w 2 = w + μ . (71)

Evaluating the integral in (71) we obtain the following cases.

Case 1: μ + C 1 = 0. Then the integration in Equation (71) gives

λ w = w + μ + C 2 ,

whence

w = 1 2 [ ( μ + C 2 ) ± ( μ + C 2 ) 2 4 λ ] . (72)

Case 2: μ + C 1 < 0. Then we evaluate the integral by writing

C 1 d w λ + C 1 w 2 = C 1 d w ( λ + C 1 ) + w 2

and write Equation (71) as

C 1 ( λ + C 1 ) arctan w ( λ + C 1 ) = w + μ + C 2 . (73)

Case 3: μ + C 1 > 0. In this case Equation (71) becomes

C 1 2 λ + C 1 ln λ + C 1 + w λ + C 1 w = w + μ + C 2 (74)

when w 2 < μ + C 1 , and

C 1 2 λ + C 1 ln w + k + C 1 w λ + C 1 = w + μ + C 2 (75)

when w 2 > μ + C 1 .

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).

5. 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 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.

Acknowledgements

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 [23] [24] [25] [26] [28] [29] [30] [31] . A concise introduction to the calculus of Lie symmetry groups can be found in [27] . Readers interested in applications of Lie groups in fluid dynamics can find a wealth of information in [23] and [32] .

Definition of one-parameter groups. Let

z ¯ i = f i ( z , a ) , i = 1 , , N , (A1)

be a one-parameter family of invertible transformations of points z = ( z 1 , , z N ) N into points z ¯ = ( z ¯ 1 , , z ¯ N ) N . Here a is a real para- meter from a neighborhood of a = 0 , and we impose the condition that Trans- formation (A1) is an identity if and only if a = 0 , i.e.,

f i ( z , 0 ) = z i , i = 1 , , N . (A2)

The set G of transformations (A1) satisfying Condition (A2) is called a (local) one-parameter group of transformations in N if the successive action of two transformations is identical to the action of a third transformation from G, i.e., if the function f = ( f 1 , , f N ) satisfies the following group property:

f i ( f ( z , a ) , b ) = f i ( z , c ) , i = 1 , , N , (A3)

where

c = φ ( a , b ) (A4)

with a smooth function φ ( a , b ) defined for sufficiently small a and b. The group parameter a in the transformation (A1) can be changed so that the function (A4) becomes c = a + b . In other words, the group property (A3) can be written, upon choosing an appropriate parameter a (called a canonical para- meter) in the form

f i ( f ( z , a ) , b ) = f i ( z , a + b ) . (A5)

Group Generator. Let G be a group of transformations (A1) satisfying the condition (A2) and the group property (A5). Expanding the functions f i ( z , a ) into Taylor series near a = 0 and keeping only the linear terms in a, one obtains the infinitesimal transformation of the group G:

z ¯ i z i + a ξ i ( z ) , (A6)

where

ξ i ( z ) = f i ( z , a ) a | a = 0 , i = 1 , , N . (A7)

The first-order linear differential operator

X = ξ i ( z ) z i (A8)

is known as the generator of the group G.

Invariants. A function J ( z ) is said to be an invariant of the group G if for each point z = ( z 1 , , z N ) N is is constant along the trajectory determined by the totality of transformed points z ¯ : J ( z ¯ ) = J ( z ) .

The function J ( z ) is an invariant of the group G with Generator (A8) if and only if

X ( J ) ξ i ( z ) J z i = 0. (A9)

Hence any one-parameter group has exactly N 1 functionally independent invariants (basis of invariants). One can take them to be the left-hand sides of N 1 first integrals J 1 ( z ) = C 1 , , J N 1 ( z ) = C N 1 of the characteristic equations for linear partial differential Equation (A9). Then any other invariant is a function of J 1 ( z ) , , J N 1 ( z ) .

Invariant equations. We say that a system of equations

F k ( z ) = 0 , k = 1 , , s (A10)

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.,

F k ( z ¯ ) = 0 , k = 1 , , s (A11)

whenever z solves Equations (A10). The group G with the generator (A8) is admitted by Equations (A10) if and only if

X ( F k ) | ( A10 ) = 0 , k = 1 , , s , (A12)

where the symbol | ( A10 ) means evaluated on the solutions of Equations (A10).

If z is a collection of independent variables x = ( x 1 , , x n ) , dependent variables u = ( u 1 , , u m ) and partial derivatives u ( 1 ) = { u i α } , u ( 2 ) = { u i j α } , , of u with respect to x up to certain order, where

u i α = u α x i , u i j α = 2 u α x i x j ,

then (A10) is a system of partial differential equations

F k ( x , u , u ( 1 ) , ) = 0 , k = 1 , , s . (A13)

Furthermore, if the transformations (A1) are obtained by the transformations of the independent and dependent variables

x ¯ = f ( x , u , a ) , u ¯ = g ( x , u , a ) (A14)

and the extension of (A14) to all derivatives u ( 1 ) , 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 [ X 1 , X 2 ] = X 1 X 2 X 2 X 1 . 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 L r and is called an r-dimensional Lie algebra. An r-dimensional Lie algebra L r 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

u α = h α ( x ) , α = 1 , , m (A15)

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 J 1 , J 2 , , .

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 jamp@scirp.org

Cite this paper
Ibragimov, R. , Ibragimov, N. , Mohazzabi, P. (2017) Approximate and Invariant Solutions of a Mathematical Model Describing a Simple One-Dimensional Blood Flow of Variable Density. Journal of Applied Mathematics and Physics, 5, 1335-1354. doi: 10.4236/jamp.2017.56111.
References
[1]   Einav, S. and Stolero, D. (1987) Pulsatile Blood Flow in an Arterial Bifurcation. Medical & Biological Engineering & Computing, 25, 12-17.
https://doi.org/10.1007/BF02442814

[2]   Elad, D. and Einav, S. (1989) Simulation of Airway Closure during Forced Vital Capacity. Annals of Biomedical Engineering, 17, 617-631.
https://doi.org/10.1007/BF02367466

[3]   Rotman, O., Shav, D., Raz, S., Zaretsky, U. and Einav, S. (2013) Biomechanical Aspects of Catheter-Related Thrombophlebitis. Journal of Biomedical Science and Engineering, 6, 6-13.
https://doi.org/10.4236/jbise.2013.612A002

[4]   Brand, M., Avrahami, I., Einav, S. and Ryvkin, M. (2014) Numerical Models of Net-Structure Stents Inserted into Arteries. Computers in Biology and Medicine, 52, 102-110.

[5]   Gummel, E., Milosevic, H., Ragulin, V., Zakharov, Y. and Zimin, Y. (2014) Motion of Viscous in homogeneous Incompressible Fluid of Variable Viscosity. Zbornik rabova konferencije MIT 2013, Beograd, 267-274.

[6]   Geidarov, N., Zakharov, Y. and Shokin, Y. (2011) Solution of the Problem of Viscous Fluid with a Given Pressure Differential. Russian Journal of Numerical Analysis and Mathematical Modelling, 26, 39-48.

[7]   Milosevic, H., Gaydarov, N. and Zakharov, Y. (2013) Model of Incompressible Viscous Fluid Flow Driven by Pressure Difference in a Given Channel. International Journal of Heat and Mass Transfer, 62, 242-246.

[8]   Dolgov, D. and Zakharov, Y. (2015) Modeling of Viscous Inhomogeneous Fluid Flow in Large Blood Vessels. Vestnik Kemerovo State University, 2, 30-35. (In Russian)

[9]   Kim, J., Ho Yoon, L. and Sehyun, S. (2015) Advances in the Measurement of Red Blood Cell Deformability: A Brief Review. Journal of Cellular Biotechnology, 1, 63-79.
https://doi.org/10.3233/JCB-15007

[10]   Bransky, A., Korin, N., Nemirovski, Y. and Dinnar, U. (2007) Correlation between Erythrocytes Deformability and Size: A Study Using a Microchannel Based Cell Analyzer. Microvascular Research, 73, 713.

[11]   Rideout, V. and Dick, D. (1967) Difference-Differential Equations for Fluid Flow in Distensible Tubes. IEEE Transactions on Biomedical Engineering, 14, 171-177.
https://doi.org/10.1109/TBME.1967.4502495

[12]   Quemada, D. (1977) Rheology of Concentrated Disperse Systems and Minimum Energy Dissipation Principle. Rheologica Acta, 16, 82-94.
https://doi.org/10.1007/BF01516932

[13]   Formaggia, L., Lamponi, D. and Quarteroni, A. (2003) One Dimensional Models for Blood Flow in Arteries. Journal of Engineering Mathematics, 47, 251-276.
https://doi.org/10.1023/B:ENGI.0000007980.01347.29

[14]   Westerhof, N., Bosman, F., Vries, C. and Noordergraaf, A. (1969) Analog Studies of the Human Systemic Arterial Tree. Journal of Biomechanics, 2, 121-143.

[15]   Formaggia, L., Gerbeau, J., Nobile, F. and Quarteroni, A. (2001) On the Coupling of 3D and 1D Navier-Stokes Equations for Flow Problems in Compliant Vessels. Computer Methods in Applied Mechanics and Engineering, 191, 561-582.

[16]   Formaggia, L., Nobile, F., Gerbeau, J. and Quarteroni, A. (2002) Numerical Treatment of Defective Boundary Conditions for the Navier-Stokes Equations. SIAM Journal on Numerical Analysis, 40, 376-401.
https://doi.org/10.1137/S003614290038296X

[17]   Streeter, V., Keitzer, W. and Bohr, D. (1963) Pulsatile Pressure and Flow through Distensible Vessels. Circulation Research, 13, 3-20.
https://doi.org/10.1161/01.RES.13.1.3

[18]   Phythoud, F., Stergiopulos, N. and Meister, J. (1995) Forward and Backward Waves in the Arterial System: Nonlinear Separation Using Riemann Invariants. Technology and Health Care, 3, 201-207.

[19]   Wittberg, L., Wyk, S., Fuchs, L., Gutmark, E., Bakeljauw, P. and Little-Gutmark, I. (2016) Effects of Aortic Irregularities on Blood Flow. Biomechanics and Modeling in Mechanobiology, 15, 345-360.
https://doi.org/10.1007/s10237-015-0692-y

[20]   Quarteroni, A., Tuveri, M. and Veneziani, A. (2000) Computational Vascular Fluid Dynamics Problems, Models and Methods. Computing and Visualization in Science, 2, 163-197.
https://doi.org/10.1007/s007910050039

[21]   Ibragimov, R.N. (2007) Oscillatory Nature and Dissipation of the Internal Waves Energy Spectrum in the Deep Ocean. The European Physical Journal Applied Physics, 40, 315-334.
https://doi.org/10.1051/epjap:2007164

[22]   Ibragimov, R.N. (2008) Generation of Internal Tides by an Oscillating Background Flow along a Corrugated Slope. Physica Scripta, 78, Article ID: 065801.
https://doi.org/10.1088/0031-8949/78/06/065801

[23]   Ovsyannikov, L. (2013) Lectures on the Theory of Group Properties of Differential Equations, Novosibirsk University Press, Novosibirsk, 1966. English Trans. Ibragimov, N.H., Higher Education Press, Beijing; World Scientific, Singapore.

[24]   Ibragimov, N. (1985) Transformation Groups in Mathematical Physics, Nauka, Moscow, 1983. English Transl.: Transformation Groups Applied to Mathematical Physics, Reidel, Dordrecht.

[25]   Olver, P. (1986) Applications of Lie Groups to Differential Equations. Springer-Verlag, New York.
https://doi.org/10.1007/978-1-4684-0274-2

[26]   Bluman, G. and Kumei, S. (1980) Symmetries and Differential Equations. Springer-Verlag, New-York.

[27]   Ibragimov, N. (2009) A Practical Course in Differential Equations and Mathematical Modeling, Higher Education Press, Beijing.
https://doi.org/10.1142/7573

[28]   Lie, S. (1891) Vorlesungen uber Differentialgleichungen Mit Bekannten Infinitesimalen Transformationen. Bearbeited und herausgegeben von Dr. G. Scheffers, B.G. Teubner, Leizig.

[29]   Cohen, A. (1911) An Introduction to the Lie Theory of One-Parameter Groups with Applications to the Solution of Differential Equations. D.C. Heath, New York.

[30]   Dickson, L. (1924) Differential Equations from the Group Standpoint. Annals of Mathematics, 25, 287.
https://doi.org/10.2307/1967773

[31]   Bluman, G. and Anco, S. (2002) Symmetry and Integration Methods for Differential Equations. Springer-Verlag, New York.

[32]   Ibragimov, N. (1990) Elementary Lie Group Analysis and Ordinary Differential Equations. John Willey & Sons, Chichester.

 
 
Top