SHOCK WAVE FORMATION IN COMPLIANT ARTERIES

. We focus on the problem of shock wave formation in a model of blood ﬂow along an elastic artery. We analyze the conditions under which this phenomenon can appear and we provide an estimation of the instant of shock formation. Numerical simulations of the model have been conducted using the Discontinuous Galerkin Finite Element Method. The results are consistent with certain phenomena observed by practitioners in patients with arteriopathies, and they could predict the possible formation of a shock wave in the aorta.

1. Introduction. We address the problem of studying shock wave formation in arteries. This was first noticed byČanić et al. in [3]. Such wave could appear if heart beats were too abrupt and then the blood could be faster than the usual blood wave. From a clinical point of view, it could be identified with the so-called pistol-shot heard, with a stethoscope, in some patients with aortic insufficiency. We have initially considered the model of blood flow in non-rigid arteries proposed by Sherwin et al. in [27,28]. We have estimated when and where the first shock wave could be formed in aorta and we have linked it with certain clinical conditions. For the interested reader in this research topic, although with different and sometimes less general models, we refer to [10,21,29].
The paper is organized as follows: In Section 2 we describe the continuity and momentum equations needed to set the model. We also add the tube law for reducing the dimension of the system, by giving a relation between blood pressure and arterial amplitude. A theoretical estimation of the location of the shock wave is described in Section 3. A numerical implementation and the results of simulations are reported in Section 4. Finally, conclusions and future research lines are outlined in Section 5. Since we are going to focus in large arteries (> 1 mm) such as the aorta, we can assume that the blood behaves as a Newtonian fluid [7,19,23]. We also assume that it is an incompressible fluid, therefore density ρ and dynamic viscosity µ are constant. For the derivation of the dynamics equations we introduce the dependent variable Q(x, t) = A(x, t)u(x, t), that will represent the volume flux at a given xsection at time t. Therefore, the three state variables of the model are A, u, and p, or equivalently, A, Q, and p. The governing equations of the model will consist of two conservation equations, and a third one responsible for modelling the artery as an elastic material.

Continuity equation.
Under the fluid dynamics continuum hypothesis, macroscopic properties such as density, pressure, and velocity are taken to be welldefined at infinitesimal volume elements. Thus, fluid properties can vary continuously from one volume element to another, and they are the resulting averaged values of the molecular properties. With this, if we take a portion [0, l] of an artery as a control volume, conservation of mass is expressed as V (t) = l 0 A(x, t)dx and, hence, the rate of mass change (or volume, if ρ is constant) is ∂V (t)/∂t = Q(0, t) − Q(l, t).
Using the definition of V (t), the fundamental theorem of Calculus, and assuming that A, u and Q are smooth enough, the resulting continuity equation is: 2.2. Momentum equation. The concept of momentum of Newtonian dynamics states that the rate of change of momentum within the control volume plus the net flux of momentum out the control volume is equal to the applied forces on the control volume. Here, we have to weight with the blood density ρ, because the flux is involved. Considering the tube as axisymmetric, this can be expressed as where f represents the friction force per unit length, and α is the Coriolis coefficient. Following [28], we assume a frictionless flow, leading to a flat profile with f = 0 and α = 1. This is a reasonable choice, since the source term is one order of magnitude smaller than the effects of non-linear advection [2]. Moreover, an inviscid flow does not generate boundary layer, and then it is physically reasonable to assume such profile [14]. Summing up, combining equations (1) and (2), we can write the system in conservative form as , and subscripts denote derivatives.
2.3. The tube law. To close the system, an additional equation is necessary, either differential or algebraic, for linking the pressure with the varying amplitude due to tube wall elasticity (known as the local tube law ). Several tube law expressions have been proposed in the literature [3,12,24,22,25,27,28,30,32,33] based upon different assumptions depending on how the flow was modelled and which theory has been used (such as shell theory [34] or plate theory [5]).
In this work, we have considered the assumption of a thin wall tube, where each section is independent of the others [27,28]. This model is based on linear elasticity, using Hooke's law for continuous media, namely σ = εE, being σ the stress, ε the strain and E the Young's modulus. The stress is a physical quantity that expresses the internal forces that neighbouring particles of a continuous material exert on each other, while the strain is the measure of the deformation of the material. The Young's modulus characterizes the stiffness due to the elasticity of the material.
Let us denote the radius of an artery by R(x, t) and R 0 (x) = R(x, 0) its initial state. Here, h 0 (x) will be used to denote the vessel-wall thickness and sectional area at the equilibrium state (p, u) = (p ref , 0), where p ref is the reference pressure. We assume a cross section of a vessel with a thin wall (h R), that is isotropic, homogeneous and incompressible, and that deforms axisymmetrically with each circular cross-section independently of the others. Making these assumptions, we can express the strain as ε = R−R0 (1−ν 2 )R0 , where ν is the Poisson's ratio, the ratio of transverse contraction strain to longitudinal extension strain in the direction of stretching force. Along with Young's modulus, ν uniquely determines the properties of a (linear) elastic material. By Young-Laplace's law [16,35], assuming that there is not external pressure, we can relate the pressure with the stress as p = h0σ πR . Combining the previous expressions, we arrive to the tube law where is the parameter determining the material properties, with A 0 (x) = A(x, 0), and p ext is the external pressure. There are no problems with the denominator in this expression, since in the materials involved ν ∈ [0, 0.5]. Moreover, although in capillaries the amplitude A 0 is very small, in those cases the wall thickness h 0 is still smaller than A 0 , and h 0 /A 0 remains bounded. on the study of the solution and its derivative along the characteristics, see [17,18]. The system (3) can be rewritten in its quasi-linear form as with To see what kind of problem we are dealing with, we look at the eigenvalues of H, that can be diagonalized as H = P −1 DP with Now, λ 1 = u + c > 0 and, since β is usually considerably bigger than the typical blood velocities, we have that λ 2 = u − c < 0. For typical values of the parameters, we refer to the compilation in [26,Appendix]. Hence, we have two real and distinct eigenvalues for the quasi-linear problem, which means that our problem is strictly hyperbolic [31]. We look for the characteristic variables (Riemann invariants) that would satisfy (5). The left and forward characteristics emanating from the origin x 1 (t), x 2 (t) can be obtained and also the characteristic variables W 1,2 , as indicated in [27,28], with Since in a physically meaningful solution β is always positive, the variables (A, u) can be written in terms of (W 1 , W 2 ).
We shall assume that U (x, t) ∈ R 2 and F : R 2 → R 2 is a smooth function of U. Since the system in (5) is strictly hyperbolic, the system in characteristic form is where W 1 , W 2 are the unknown functions and λ 1 , λ 2 are smooth functions of W 1 and W 2 . Due to the hyperbolicity, this can be done at least locally. We also assume that the system is non-linear in the considered domain, that is ∂λ1 Here we have the initial boundary-value problem where we can assume without loss of generality that x 1 (0) = 0. Canić and Kim also found sufficient conditions for the hyperbolicity of the system [2, Th. 3.1]. They decomposed the domain as Figure 2.
They first showed the hyperbolicity on the bounded domain D T 2 , T > 0, which is the restriction of D 2 to the semiplane of t ≤ T . Later, sufficient conditions were provided for the existence of a unique C 1 solution (W 1 (x, t), W 2 (x, t)) of the initial boundary problem given by (10)-(11) on D. However, in our case one of those conditions does not hold, since we have λ 2 = W 2 , then ∂λ2 ∂W2 = 1 > 0.  The question that arises is whether the lack of smoothness could be plausible. In our main result we obtain an estimation of the first shock wave appearance. We follow the steps ofČanić and Kim [2, th. 3.1], applied to the model by Sherwin et al. [28]. Keener and Sneyd also obtained a similar (but less general) result in [15]. For the sake of clarity we use the notation u t = ∂u ∂t in our development. Theorem 3.1. Let us consider the initial boundary problem stated in (10)- (11). Let us assume constant initial data A(x, 0) = A 0 , u(x, 0) = 0. Let ω be the first time when the forward characteristic intersects the left spatial boundary. Then, the the time t s of the first shock formation is given by Proof. In terms of Riemann invariants, the initial data read W 0 1,2 (x) = ±4 β 2ρ A 1/4 0 . Then, W 1 is constant everywhere in the region of smooth flow D 1 . The characteristics x 2 (t) are straight lines in D 2 emanating from (0, 0), with x 2 = λ 1 . Therefore, the solution on D 1 is bounded by the left boundary x 1 = 0 and the forward characteristic x 2 , see Figure 2.
There, the solution is driven by u(·, t) on x 1 (t) and will develop shock waves due to the fact that u t (·, t) changes of sign. We use the following expressions of u(x, t) and A(ω, t) taken from [2, Eq. 3.15 & 3.16], in which x (ω, t) is the forward characteristic passing through the point (ξ, ω). with To estimate the time and location of the shock wave formation (t s , x s ) we take into account that the partial derivative ∂W 1 /∂x at (t s , x s ) must blow up. This occurs at the point where the denominator in (13) vanishes. Hence, the time t s can be calculated by recalling that λ 1 = W 1 , which implies that ∂λ 1 /∂W 1 = 1, and that W 1 = W 0 1 everywhere. This implies in (14) e h(W1(ξ,ω),W2(ξ,ω))−h (W1( x(ω,τ ),τ ),W2(ξ,ω)) = 1 (15) and then Differentiating and then we can conclude Since ∂W 2 /∂x = 0, we obtain Therefore, isolating, we get the first time t s in which a shock wave is formed An inspection of (20) indicates that the shock will be produced sooner if the inflow u(0, t) accelerates or/and if the walls of the vessel are less rigid (due to the β factor). In order to explore the connections of the shock wave with the pistol-shot heard in aortic insufficiency, we take as an approximation the measures in [11,19,6]. We assume A 0 = A(0, t) ≈ 4×10 −2 m and the blood flow velocity as 1 m/s, according to [20]. Also β is computed with the parameters in [26,Appendix]. For a healthy person we have taken u t (0, t) = 7 m/s 2 , following the correlations of [20], and a Young's modulus of E = 10 5 Pa. Hence, the values of time and place for the first shock wave are t s ≈ 1 + 4 29633.5/(2 × 1050)(4 × 10 −2 ) 1/4 7 ≈ 1.1s, x s = t s λ 1 ≈ 8.5m, (21) a distance that, according to [9], is far from the mean length of aorta, 33.2 cm. Now, for a patient with aortic insufficiency, the heart increases its volume and, since the aortic valve does not close properly, the muscle must do a greater contraction. As a consequence, a greater blood flow acceleration happens in each heart beat, see [13]. Using the aforementioned bibliography, we take a value of u t (0, t) = 15 m/s 2 . Moreover, in some cases aortic insufficiency comes associated to a decrease in wall rigidity (also associated to aneurysm), which could be represented by taking a value of E = 2 · 10 3 Pa. With this new choice of parameters, we have: 13s, x s = t s λ 1 ≈ 0.25m. (22) This result indicates that the formation of shock waves inside the aorta is consistent with the model. Next, we present some numerical simulations that describe our progress in the simulation of shock waves. 4. Numerical simulations. In this section we provide the results of a numerical simulation using a set of parameters compatible with the formation of a shock wave. The propagation of the velocity and the amplitude have been simulated using the Discontinuous Galerkin Finite Element Method (DG-FEM) as in [27]. For being applied, a modal discretization of the solution has been done. Here, the state variables are not values of the solution, but its coefficients in a certain basis. In our case, we have considered the base of Legendre polynomials. Therefore, at each time step we compute the coefficients of each polynomial up to certain degree. The numerical solution given by a DG-FEM scheme can be discontinuous at element interfaces and this discontinuity is resolved by the use of a so-called numerical flux function, which is a common feature within Finite Volume Methods. We refer to [8,26] for a more detailed explanation of the implementation. This method is proposed by Sherwin et al. [28] for its numerical robustness when a smooth flow is present. However, depending on the degree used for polynomials, it can provide a solution too smooth to reproduce strong shock waves. As a consequence, we only expect to be able to observe weak waves.
In their simulations, Sherwin et al. [28] show the feasibility and usefulness of the method for this type of problems using some non-physiological parameters to augment the effect of the travelling wave. In our experiments we have chosen a set of parameters that are more physiologically meaningful.
Next, we describe the experimental setup. Let T be the period of a wave and w = 2π/T . We consider a beat-like left inflow velocity as with 1 [0,T /2] (t) equal to 1 if 0 ≤ t ≤ T /2, and 0 otherwise, see Figure 3. In the present simulation, we have chosen the stiffness parameter β = 10 4 , the period T = 0.8, a blood density of ρ = 1050 kg/m 3 , both null initial blood flow velocity and external pressure, and an inflow artery amplitude of 4 · 10 −2 m. The inflow and outflow amplitudes coincide with the initial one. The outflow velocity also coincides with the null initial one. Although this assumption could affect the simulation when the traveling wave reaches the end of the artery, our simulation only reproduces the initial phase of the wave formation when it travels along the first elements, far from the outflow boundary condition.
For the DG-FEM discretization we use a cubic polynomial basis in the computational interval of [0, 200], with 10 equally distributed elements. The time step is ∆t = 10 −3 . In Figure 4 we can appreciate the progressive formation of the weak shock wave in the first element of the mesh. Note that he discontinuity present at x = 20 corresponds to the frontier between two elements and not to a discontinuity in the solution. This discontinuity appears as a result of the DG-FEM formulation, that provides two values for the velocity in inter-element frontiers. We note how this does not appear as an abrupt anomaly in the velocity function, but as a slight "accumulation" when the weak shock is about to be formed. We can appreciate this weak shock more clearly in the last plot of figure 4 (146 ms) around 18-19 cm. This weak shock should not be seen as a classical shock wave due to its magnitude, but as a zone of accumulation of the velocity. This also makes clinical sense since if these shock waves were too sharp and sudden they will cause further complications resulting in more morbid conditions. In Figure 4 we can also see how the magnitude of the inflow velocity (rising in the top right plot) allows the shock wave to be formed, although not instantly. In a smoother framework we could observe how the transition between times 80 and 146 conserves the crest of the wave. Based on these numerical results we can say that the traveling wave evolves to a weak shock wave in the range of time and distance simulated. It has to be still analyzed whether we will be able to simulate strong shock waves with the numerical method used or if, on the contrary, the smoothness of the numerical solution makes it difficult to show it without using very high degree polynomials. In the later case, we consider using methods specific for shock wave detection, such as schemes used in similar studies on nonlinear acoustics [4]. 5. Conclusion. We have started the study of finite time shock wave formation in arteries under the compliant artery 1D-model by Sherwin et al., showing that, after a time interval given by (12), a shock wave appears. This extends previous work in shock wave formation byČanić et al. [3].
We have considered a set of parameters closer to the physiological interpretation than in previous works on the topic. Our results predict the appearance of a shock wave inside the aorta for parameters compatible with the existence of arterial insufficiency. Nevertheless, in order to confirm the connection with the pistol-shot, a deeper and more precise study of the parameters should be carried out.
The model has been implemented using Discontinuous Galerkin FEM, and this has permitted us to reproduce the appearance of a weak shock wave while preserving the smoothness of the flow. An interesting point to be addressed is to determine what is the range of parameters that still can be used for simulations with this method and whether other numerical methods will improve the simulation of the shock wave [4]. Other possible research line is to study the properties of blood flow in compliant arteries compared to flow of other types of fluids in deformable pipes [1].