The 2-component dispersionless Burgers equation arising in the modelling of blood flow

This article investigates the properties of the solutions of the dispersionless two-component Burgers (B2) equation, derived as a model for blood-flow in arteries with elastic walls. The phenomenon of wave breaking is investigated as well as applications of the model to clinical conditions.


Introduction
There are many examples of systems of nonlinear PDEs modelling the propagation of waves in fluids under various conditions. One of the earliest models was of course the KdV equation, modelling the propagation of surface waves in shallow water. Meanwhile the nonlinear Schrödinger equation has proved very successful at modelling the propagation of wave fronts in deep water.
The Camassa-Holm equation [7] u t − u xxt + 2ωu has gained popularity as an integrable model with many applications and interesting properties [5,11,6]. Among its many applications, it is well known as a model of shallow water waves, admitting many novel solutions cf. [27,14]. Another physical application of the system arises in the study of axially symmetric deformation waves propagating in hyperelastic rods cf. [17,16]. The system also admits various generalisations, the most popular of which is the two-component one (CH2) [32]: where m = u−u xx . Taking ρ = 0 reduces the CH2 system to the CH equation (1). The CH2 energy Hamiltonian is given by and is clearly positive-definite. The CH2 system (2)-(3) is bi-Hamiltonian. This means it possesses two compatible Poisson brackets. The first Poisson bracket between two functionals F and G of the variables m and ρ is in semidirect-product Lie-Poisson form [32,12,26]: This Poisson bracket generates the CH2 system from the Hamiltonian H 1 . Its second Poisson bracket has constant coefficients, and corresponds to the Hamiltonian H 2 = 1 2 (uρ 2 + u 3 + uu 2 x )dx. There are two Casimirs for the second bracket: ρ dx and m dx.
The CH2 model has various applications. For example, in the context of shallow water waves propagating over a flat bottom, u can be interpreted as the horizontal fluid velocity and ρ is the water elevation in the first approximation [12,26]. In the shallow water regime there are scale characteristics, one of which is δ = h/λ where h is the average depth and λ is the wavelength of the water waves. With these scale factors, if u = O(1), then m = u−δ 2 u xx , see e.g. [12,26]. In the limit λ >> h or δ → 0, we have m = u, which we call the 'two component dispersionless Burger's equation': The first Poisson bracket (5) remains valid in the case m = u with Hamiltonian H 1 . The second Poisson bracket (6) in this limit becomes and the corresponding Hamiltonian is H 2 = 1 2 (uρ 2 + u 3 )dx. Dispersionless systems like (7)-(8) also arise in the context of very long water waves, in particular in the modelling of tsunamis as they approach the shore [13].
Both systems are illustrated in Fig.1 and Fig.2 with the so-called 'dambreak' initial conditions. The dam-break problem involves a body of water of uniform depth, initially retained behind a barrier. When the barrier is suddenly removed at t = 0, the water flows downward and outward under gravity. The problem is to find the subsequent flow and determine the shape of the free surface. This question is addressed in the context of shallow-water theory, e.g., by Acheson [1], and thus serves as a typical hydrodynamic problem of relevance for shallow water equations. Dam break initial conditions are: solved in a periodic domain of width 100 for the CH2 and Burger's equations. Both are solved in conservative form using an implicit midpoint rule. For CH2 the momentum equation is written as For Burger's equation, a potential form is used In Section 2 we develop a model of blood flow in elastic tubes modelling arteries. The blood is modelled as a Newtonian fluid, which is a good approximation in the case of the main arteries, as the diameter of the artery is much larger than the cross-section of a typical blood cell. In Section 3 we show that the model developed in Section 2 may be solved directly by using characteristic curves. Section 4 examines a related model, however unlike the primary model we work with, the model examined in Section 4 does not possess a positive definite Hamiltonian. We find explicit solutions to this system based on a method analogous to the method o characteristic outlined in Section 3. In Section 5 we show that the two component model developed in Section 2 does display the phenomenon of wave breaking, and we give a physiological interpretation of the result based on clinical observations. In this section we also discuss the initial conditions for which the system possess global solutions.

Modeling Newtonian fluids in elastic tubes
The 2-component Burgers equation, and related systems, are well understood quasi-linear systems arising in many physical applications, for example as model of shock-waves in gas dynamics (see e.g. [34], Section 13.2). In this paper we give a derivation of the system (7)- (8) in the context of a Newtonian fluid within an elastic tube, modelling the flow of blood within arteries.
The model is quasi-linear in the dynamical variables, u(x, t)-the fluid velocity, and A(x, t)-the tubes cross section. The physiological relevance of the current article is in relation to the so called pistol-shot pulse. This is a popular term for a phenomenon reported by clinicians in which a loud cracking sound is heard by the stethoscope over an artery, caused by a large distension followed by an abrupt collapse of the artery wall. This phenomenon is known to occur in large arteries during aortic regurgitation, a condition in which some blood leaks back into the left ventricle during systole.
To establish a model of arterial blood flow, we consider a simplified model of blood itself, in that it will be modelled as an incompressible, inviscid and irrotational Newtonian fluid. The artery is modeled as an axially symmetric elastic tube, with cross section A(x, t) which depends on the axial coordinate x and time t. Being elastic, the artery experiences a restoring force when expanded from it's equilibrium cross-section, and as such exerts additional pressure on the blood contained within.
The volume of blood contained in the infinitesimal volume of artery between the axial locations x and x + dx, at a fixed instant t is dV = A(x, t)dx. In addition, the mass content in that same volume of blood will be the material density at the location at that instantρ(x, t), times the volume itself, or The sub-script (t) on M (t) denotes the mass of the blood at a fixed instant t.
Since we model the blood as an incompressible fluid, the density of the blood is homogenous throughout the artery at all instants, and so in our model the blood density is a fixed constantρ. Within a finite volume of artery, between the axial locations a and x, with a < x, and at some fixed moment t, it follows that the mass content is, The sub-script appearing in V ax denotes the volume between the axial locations a and x.
The velocity of the blood through the artery at a fixed axial coordinate and time will be uniform throughout the cross-section at that location and time. In an infinitesimal time element dt, the blood content at the location x will be displaced by a distance dX = u(x, t)dt. In this time interval, the total mass of blood to traverse the cross-section A(x, t) will be It follows that the change in the blood content in the arterial volume between a and x, in the time interval dt, is the total blood displacement into the volume across A(a, t) less the total blood displacement out of the volume across A(x, t). So we may write, No subscript appears on the mass element dM on the left hand side above, since we are evaluating the change in mass content over the infinitesimal time interval dt.
We may rewrite the above equation as one involving the rate of change of blood content within a finite arterial volume as follows, Operating on this equation with ∂ x , the fundamental theorem of calculus applied to the left-hand side gives us while ∂ x applied to the right hand side gives us Thus we arrive at the continuity equation, where sub-scripts denote partial differentiation with respect to the particular variable. The factor ofρ, common among all three terms has been cancelled from the above equation.
The physical interpretation of this equation is as follows. In an infinitesimal time interval dt, the infinitesimal increase (decrease) in the blood content dM at x causes an increase (decrease) of the blood volume, since the blood itself is incompressible. As the walls of the artery are elastic, this increase (decrease) in blood volume is accommodated by an expansion (contraction) of the artery, which is realised as an increase (decrease) in its cross section A(x, t), at the location x.
Newton's laws dictate that any forces acting upon the blood will cause a corresponding change of momentum of the blood. In the case of fluids, the forces acting upon any individual element may be separated into two categories • The external forces acting on the fluid elements, which in this case will be the restoring forces present in the artery walls, thereby exerting an external pressure on the blood at a given location x and some instant t.
• The internal forces acting upon a fluid element caused by contiguous fluid elements.
Both external and internal forces are found to cause an overall change in the linear momentum of an individual fluid element, in accordance with Newton's second law. The linear momentum of an infinitesimal volume of blood at some location x and instant t isρu(x, t)dV x =ρu(x, t)A(x, t)dx, or the infinitesimal mass element times the velocity thereof. The time rate of change of the momentum, namelyρ∂ t (u(x, t)A(x, t))dx, will appear as the left hand side of Newton's second law. In the case of a finite volume of blood this becomes which is the rate of change in momentum of the blood contained within the artery between a and x at the instant t. As we have already seen, an infinitesimal time displacement causes a corresponding spatial displacement of the blood at x, by an amount dX. However, in moving this distance the velocity of that blood will change by the amount It follows that the corresponding change in linear momentum of this blood will be it mass times the change in velocity, or which for a finite volume becomes Thus we have established the change in linear momentum of a finite volume of blood, due to the spatial displacement thereof, during an infinitesimal time increment dt.
Next we must also include the possibility that the linear momentum of the blood will change in an infinitesimal time dt, due to a change in the blood-content. Suppose at some instant t the rate of increase of blood at x is u(x, t)A(x, t), while the rate of decrease of blood at x + dx is It follows that the rate of change of the mass will be −ρ(u(x, t)A(x, t)) x dx, and the corresponding rate of change of the linear momentum will be −ρ(u(x, t)A(x, t)) x u(x, t)dx.
In the case of a finite volume, the corresponding change in linear momentum is, This concludes the effect of internal forces acting upon an individual fluid element within the artery.
Next we need to investigate the effects of external forces on a fluid element. In this case the external forces are provided by the restoring forces within the arterial wall itself, which exerts a force on the blood. Between the axial locations x and x + dx, the internal surface of the artery has a directed area elementndS, withn being the outward normal and dS being the magnitude of the infinitesimal area element. The force exerted by the tube on the blood is −P (x, t)ndS, where P (x, t) denotes the pressure exerted due to the restoring force of the artery acting on the blood. In the case of a finite volume V between a and x, we find the corresponding force exerted by the artery with surface area S to be − S P (ξ, t)ndS.
Applying Gauss' law this may be re-written as This accounts for all the forces acting on the blood contained within the artery between axial locations a and x.
Having established the rate of change of momentum, along with the various contributing forces, we require an equation of motion for this finite fluid element. On applying Newton's second law to the above acceleration and forces, we obtain, Dividing both sides byρ, operating on the resulting equation with ∂ x , and applying the fundamental theorem of calculus to both sides we obtain, Expanding the derivatives on each side, applying the continuity equation (12), and dividing by A(x, t), we find We see that in our model the blood flow will satisfy Euler's Equation for an incompressible, inviscid and irrotational fluid.
In general it is difficult to obtain from first principles an explicit expression for the transmural pressure P (x, t). However a large body of experimental data exists to suggest a plausible correspondence between the pressure and the cross section of the tube itself. Indeed, a standard example is the so called Windkessel model [28], in which the pressure is related linearly to the cross-section, The constant k W relates the elastic restoring force of the tube when distended to cross-section A(x, t), to the pressure exerted on the blood. This is generally determined from clinical data. In this paper we will adopt a slightly more complicated model, whereby the relationship between pressure and aortic cross-section is quadratic, namely P (x, t) = P (A(x, t)) = kA 2 (x, t).
Again, the constant k in our model is determined from clinical data, however it shall play no significant role in any of the remaining results. As such, we shall redefine our cross-section, such that k ρ A 2 (x, t) → 1 2 A 2 (x, t). Finally the second governing equation of our model for arterial blood flow is, The equations (12) and (14) when taken together constitute a two-component dispersionless Burgers equation.

Solutions when k > 0
In this section we consider in more detail the behaviour of the system, where we have chosen a definite sign for the physical parameter k = 1. As we already mentioned, one of the great advantages of this system is that it may be solved directly by the method of characteristics. To illustrate this we shall define a pair of diffeomorphisms ψ ± (x, t) by the following criteria: When we apply ∂ t once more to (16), and using (15) we find, Since ∂ 2 t ψ ± (x, t) = 0, it follows that ∂ t ψ ± (x, t) must be constant in time, and may depend on x only.
Using this, and the definition supplied in (16), we see that the physical variables evaluated along the flow of ψ ± (x, t) must satisfy the constraints, We see that (17) requires ψ ± (x, t) be linear in t, while the second equation in (16) imply that ψ ± (x, t) satisfy, Applying ∂ t to ψ ± (x, t) as it is given in this expression, and comparing to (16) and (18), we can see that the functions γ (±) 0 (x) may be written in terms of the initial data, γ (±) 0 (x) = u 0 (x) ± A 0 (x). So we see that we may solve (16) to find, allowing us to express ψ ± (x, t) in terms of the given initial data u 0 (x) and A 0 (x). Since are diffeomorphisms, at least for appropriately chosen initial data, it is in principle possible to invert (19), and upon doing so one may obtain explicit solutions for the physical variables u(x, t) and A(x, t) in terms of the initial data (u 0 (·), A 0 (·)) : where it is understood that ψ −1 ± := ψ −1 ± (x, t).
Example: As an example, we will consider a solution in which our initial data is of the form u 0 (x) ∼ A 0 (x) ∼ x 1 3 . Specifically, our initial data is defined by, where a ± are constants. It follows from (19) that the diffeomorphisms ψ ± (x, t) may be written as, Making the invertible substitution x → w 3 , and with the corresponding change of variables ψ ± (x, t) → y ± (w, t), we may re-write the diffeomorphisms in (22), We now have a pair of monic polynomials with argument w ∈ R, namely, the discriminants of which are given by, The values of the discriminant determine the quantity and nature of solutions for w. In particular, we are interested in real solutions w(y ± , t).
Equivalently, this solution allows us to solve for x in terms of ψ ± (x, t) and t, that is, it offers us an expression for ψ −1 ± (x, t). Depending on the discriminant D, we may have several real roos, all of which may be distinct, or some of which may be equal. In the case D = 0, all roots are equal and real, for a ± ∈ R and real ψ ± (x, t). Moreover, the discriminant allows us to explicitly calculate the time at which the wave breaking occurs, In this case, the functions ψ ± (x, t) do not have unique inverses, and so no longer behave as diffeomorphisms. Furthermore, the corresponding solutions u(x, t) and A(x, t),as defined by (20), will no longer satisfy |u x (x, t)| < ∞ and |A x (x, t)| < ∞, since ψ ± (x, t) are no longer strictly monotone increasing functions of x, for all t > 0. Remark The functions ψ ± (x, t) are diffeomorphisms only in the case where wave-breaking does not occur. We assume (u 0 (x), A 0 (x)) ∈ C 1 × C 1 and also, u 0 (x) and A 0 (x) are bounded, that is, Therefore, our solutions blow-up only if u ′ 0 (x) ± A ′ 0 (x) < 0 at some point and t > 0, otherwise the solutions are global.

Solutions when k < 0
A related system is the so called wrong sign Burgers equation given by, which differs from (7,8) by the sign of the (A 2 ) x term. In this case we set k = −1. It can be shown that the Hamiltonian for the system may be written as, which we can see does not have a definite sign. It follows that the system appearing in (26) will describe physically unstable systems. In analogy to the case of solution via characteristics, we may construct a pair complex conjugate mappings χ ± (x, t), which are formally defined by, As in the case of (18), we find the analogous relations, The construction of such a pair of solutions was given in [29], wherein the authors obtained solutions by the method of Riemann invariants, with the mappings χ ± (x, t) given by It follows that the initial data corresponding to each of these mappings is given by, with the parameter λ being freely adjustable. The corresponding solutions are found to be, where the solution u(·, t) is the real root of the cubic polynomial, for arbitrary ξ ∈ R.
Having an explicit expression for A(χ ± , t) in term of u(χ ± , t) we now require an explicit expression for u(χ ± , t) in terms of x and t. To proceed, we notice from (29) that our solution u(χ ± (x, t), t) must satisfy It follows from (32) and (33) that the solutions u(χ ± (x, t), t) may be written as, Substituting this expression into the cubic polynomial in (34), we find, and so we also have explicit solutions for A(χ ± , t) in terms of x and t, as follows from (33).

Solutions that are not continuous diffeomorphisms
The phenomenon of wave-breaking was one of the most interesting and exciting aspects of the Camassa-Holm equation to be investigated following its discovery in the work [7]. In this article we have focussed on a similar though somewhat simplified system to model the flow of blood in arteries, and we would like to investigate the possibility of such a phenomenon arising in the context of this model. The clinical relevance of this is in relation to the phenomenon of the pistol shot pulse. Such behaviour arises during systole of patients with aortic insufficiency, when blood ejected into the artery is regurgitated back through the aortic valve, into the ventricles. To maintain systole pressure, the heart will contract more during ventricular systole, thereby ejecting blood into the artery with greater pressure. After systole, the ejected blood will induce a pressure gradient along the radial artery, however, owing to the initial aortic regurgitation, there will be insufficient blood mass to maintain the pressure throughout the length of the artery. This in turn will cause a sudden increase followed by a rapid decrease in the aortic cross section, which is observed as the femoral pistol-shot pulse.
In this section we aim to establish the conditions under which a sudden expansion in the aortic cross section is followed by a collapse thereof, in rapid succession. We aim to show that the phenomena of wave-breaking arising in the system (7)- (8), which we are using as a simplified model of aortic blood flow, is sufficient to account for this clinical phenomenon. Mathematically speaking, wave breaking occurs when our solutions u(x, t) and A(x, t) remain bounded for all (x, t) ∈ R × [0, T ), while the magnitude of their gradients become singular in finite time [10,35], We begin by constraining our initial data, such that, We will show that continuous solutions u(x, t) and A(x, t) do indeed remain finite, if their slopes remain finite. To do so, we first multiply each member of the system in (15) by u(x, t) and A(x, t) respectively. Then integrating over R and imposing the boundary conditions, we find the following conditions Indeed we see from this result that the quantity, is an integral of motion, and may act as a Hamiltonian for our system. Next we operate with ∂ x on each member of (2), and multiply by u x (x, t) and A x (x, t) respectively. Imposing the boundary conditions (38) and integrating over R, we find, Using (39) and (40), we find that, It follows from (41) that where sup Gronwall's inequality [4] states that for a pair of functions f (t) and g(t), continuous on some interval α ≤ t < β, and with f (t) differentiable on (α, β), and such that Upon applying Gronwall's inequality to (42), we find that Here we introduce 0). In addition, for any function u(x, t) ∈ H 1 , we have |u(x, t)| 2 ≤ u 1 = R (u 2 (ξ, t) + u 2 ξ (ξ, t))dξ.
Since we assumed (u(x, t), A(x, t)) ∈ H 1 (R) × H 1 (R), it follows from this inequality, along with the result in (43) that It follows that the solutions |u(x, t)| and |A(x, t)| remain bounded for all (x, t) ∈ R × [0, T ) if |u x (x, t)| and |A x (x, t)| remain bounded for the same values of x and t. Next, we would like to establish under what conditions the functions u x (x, t) and A x (x, t) actually become singular, while u(x, t) and A(x, t) remain bounded. We return to the diffeomorphisms introduced in Section 3, in particular (18). Differentiating this once with respect to x, we find Substituting (19) into this, we find .

Conclusion
We have derived a model for the flow of blood through an artery. We model the blood as a Newtonian fluid, flowing through an elastic tube. The restoring forces in the elastic tube are assumed to depend on the distension or contraction of the artery from its equilibrium cross section. Specifically, the pressure exerted on the blood by the artery is proportional to the square of the cross-section.
One particular nonlinear effect we investigate in the model is the steepening of the wave front as it moves away from the ventricle during systole. The same phenomenon was studied in other models by several authors, in particular in the works [2,3,30].
In the event that the wave front becomes too steep, the top of the wave front overtakes the bottom, and a shock (or discontinuity) develops as a solution, which is typical of hyperbolic equations. Of course a true shock is not possible in this context as blood viscosity and the elastic properties of the arterial wall preclude the formation of a discontinuous solution. Nevertheless, it is possible to generate very steep pressure gradients within the artery, due to aortic regurgitation during systole. Under such circumstances, clinicians report the phenomenon of the pistol shot pulse, in which the arterial wall undergoes an abnormal expansion, followed by a sudden collapse, which we interpreted here as the development of shock in the wave-front [28].