Numerical Solution for the Boussinesq Type Models with Application to Arterial Flow

A family of classical Boussinesq system of nonlinear wave theory is presented in a form of conservative equations supplemented by Riemann solver, which is a fundamental block in the Godunov frame formulation of flow problems. The governing system simulates different physical phenomena, such as propagation of a small amplitude waves on a surface of water, or pulse wave propagation in compliant thin walled arterial systems. While the first model is for verification purpose only, pulse wave propagation through the junction of thin walled elastic branches is of primary interest. The inertial effects associated with transversal motion of the wall are introduced, affirming the dispersive nature of pulsating waves. The problem of accounting for branching and possible discontinuity of wall properties is addressed. Preliminary analysis is presented which leads to the correct jump conditions across bifurcated area. As a result, the model accurately describes the formation of transmission and reflection waves at bifurcation, including effects of discontinuities. The Riemann solver supplies the inter cell flux and junction fluxes, based on conservation of volume, momentum and energy, with account of losses associated with the flow turn angle at bifurcation. An implicit monotonic total variation diminishing (TVD) scheme, second order accurate in time and space has been applied for the analysis of solitary wave solution in bifurcated arteries. Numerical results are in a good agreement with the known analytical and numerical solutions reported elsewhere. Based on direct computational analysis the inverse solution was obtained, calculating the local elastic properties of the arterial wall, using typical diagnostic measurements. Mathematical modelling presented in this work leads to a physiological understanding and interpretation of diagnostic measurements of the wave forms of a blood pressure, flow rate and an artery wall deflection. Unrestricted use, distribution, and reproduction in any medium are permitted, provided the original work is properly cited. Nomenclature TVD Total variation diminishing CFD Computational fluid dynamics PWV Pressure wave velocity c speed of propagation cMK Moens-Korteweg speed of propagation A Cross sectional area í µí¼ Cross sectional axial flow velocity í µí¼ Circumferential wall strain p Static pressure í µí¼ Fluid density E Elastic modulus h Wall thickness r radius Au Peak flow velocity AQ Peak volume flow rate Ap Peak pressure


Introduction
Consider a broad family of Boussinesq type dispersive wave models [1] presented in a conservative form where dot above variable means time derivative, subscripts are partial derivative by the axial coordinate.The equations (1,2) describe a diverse range of physical phenomena, for example acoustics models [2], nonlinear Korteweg de Vries or Benjamin-Bona-Mahony dispersive waves models [2,3], propagation of small amplitude waves on the surface of incompressible flow in the passage of the finite depth [2], pulse propagating incompressible flow in elastic or viscoelastic segments [4][5][6][7], and others [2,3].
A well-developed high resolution computational fluid dynamic (CFD) technique, capturing the physics of a wave structure by applying the relating Riemann solver, is presented.Convective terms have been discretized by the third order implicit finite volume approximation, whereas diffusive and dispersive terms have been approximated based on second order finite difference implicit scheme.To avoid numerical oscillations across discontinuities of solution, a total variation diminishing (TVD) scheme was applied [5].The proposed employment of implicit solver is justified by significant restriction imposed on a time step by alternative explicit scheme in case of stiff boundary value problems.
The problem of accounting for topological bifurcation and possible discontinuity of wall properties is addressed.Preliminary analysis is presented which leads to the correct jump conditions across bifurcated area.The Riemann solver supplies the inter cell flux and junction fluxes, in case of bifurcation, upheld by conservation of volume, momentum and energy flow rates.Based on direct computational analysis the inverse solution was obtained, calculating the local elastic properties of the arterial wall using typical diagnostic measurements.

Dynamics of Incompressible Flow in Elastic Vessel
Conservation of mass and momentum results in the following system of one dimensional equations (A=A(x,t)-cross sectional area of the vessel, u=u(x,t)averaged flow velocity, p=p(x,t) -static pressure, f -the source friction term, ρ -density of incompressible fluid) [4].
For the thin walled elastic membrane vessel in a transversal motion, pressure -strain relationship is maintained by dynamic equilibrium condition ( w  - density of the wall, , where E elastic modulus, ν-Poisson coefficient, η circumferential strain, h thickness of the vessel, r0 radius of unstressed vessel) , the total system of equations can be presented in the following nonconservative form

Solitary Waves Solution
It is well known that Boussinesq equations possess solitary wave solution [1,2].Substituting solution in a form of a traveling wave (, ) = (), (, ) = (), (, ) = (), where =  −  , into (6),( 7), (8), and integrate, one can obtain ( , where prime means derivative by X) A non-linear momentum equation is obtained by multiplying equation (10) by  d and integration to yield Setting Aη=maxX(η(X)), AP=maxX(p(X)), Au=maxX(u(X)), and using (9) to calculate integrals in equation (11), one can find speed of a wave propagation (pressure wave velocity (PWV)) in a form (12) where: 12) presents the pressure wave velocity of the forward and backward solitary wave propagation.Based on asymptotic expansion , it can be reduced to the simplified form The increment  9), (10) for one way wave propagation   =  − ( −  0 )/(1 +   ) 2  (15) shows dependence of wave velocity propagation on wall deflection and pressure amplitudes.All quantities are presented in a non-dimensional form, being normalized to cMK for all velocities, to  ×   2 for the pressure, to the radius r0 of the unstressed tube for the normal deflection.
Results show that the speed of propagation strongly depends on the amplitudes of pressure waves.Not mentioning theoretical importance, equations (15-17) are useful for verifying numerical accuracy of CFD based simulation of solitary wave.
Figure 1a.The dependence of non-dimensional wave velocity () propagation on peak normal deflection and flow velocity amplitudes ( 0 ). Figure 2b The dependence of non-dimensional wave velocity () propagation on peak pressure and flow velocity amplitudes ( 0 )

Numerical Approach
By inserting (5) into (3-4) we obtain the system written in the matrix conservation form The unknown vector of primitive variables should satisfy partial differential equation (18) supplemented with relating initial and boundary conditions.To develop implicit procedure of integration with respect to time we invoke an asymptotic differential equation relating the unknown vector Q to its increment with respect to time Q  on different time layers.
Here, n is the time layer number,  is the time step,  -scheme parameter.Linearizing (18), we obtain We substitute expressions for  Q and  Q given by (18), ( 19), (21), ( 22) respectively into (20) to obtain a differential equation for q  , in which the coefficients and the right hand side are calculated in terms of variables on the preceding time layer.
To construct an implicit finite difference scheme, we express the flux Jacobian matrix as the difference between two nonnegative definite matrices defined as To approximate the differential equation in computational domain, we construct a uniform grid with a grid size h.Integer index corresponds to the cell center, fractional (half integer) to the cell face.Right and left finite difference operators are defined as follows: . As a result equation ( 23) is replaced by the following finite difference counterpart The finite difference approximation of RHSi is defined as where fluxes are evaluating according to TVD procedure using limited extrapolations of primitive variables described with details by [6][7][8][9] The general Riemann problem for discontinuity breakdown is then utilized to evaluate the inviscid fluxes F on the sides of each cell.The "viscous" fluxes are calculated using standard central differences [7].The convective part of implicit operator is achieved by using the upstream differences with consideration of the sign of the characteristic speed (eigenvalues of the Jacobean matrix).Implicit boundary conditions for variations of variables are incorporated.The basic solver strategy employs a Gauss-Seidel relaxation procedure.In addition, Newton sub iterations are used to achieve convergence at each global time step reducing linearization errors.The procedure is of second order accurate in both space and time and provides stable convergence for time dependent hyperbolic equations without any artificial dissipation.Further details of numerical technique can be found in [8][9][10][11][12][13].

Numerical Experiments 3.1. Scalar nonlinear waves
Consider the one dimensional non-steady viscous Burgers' equation The Riemann problem with an initial discontinuity at x=0 can be expressed as a solution of the corresponding "inviscid" Burgers equation ( where the subscripts L and R denote the left and right regions respectively.As it is known the corresponding solution U consists of shock and rarefaction waves [2], propagating from discontinuity For the first example the Burgers' equation was solved for a few viscosity coefficients ε=0.25, 0.1, 0.05, 0.01, using 100 cells, with initial conditions , where . In Fig. 2 results are compared with exact solution found in [2].Courant number up to 30 does not practically affect an accuracy of results.A snap shot of evolving sinusoidal profile is shown in Fig. 2, where viscosity coefficient is low ε=0.003(stiff problem).As expected no oscillations are observed using TVD formulations.

Flow in junction
Consider junction element of bifurcation tree shown in Fig. 4. We denote the parent branch by an index 1, and two daughter branches by the indices 1 and 2. Variables related to the cell faces are noted by capital letters, the ones related to the cell interior -by low case letters.We have six unknown variables, pressure and velocity values associated with cell faces: and three equations describing the transfer of forward and backward propagating Riemann invariants [4] Analytical solution can be obtained based on a linearized form of the governing system (30), (31), which was used in the present work to calculate fluxes.For general purpose the linearized solution could be used as an initial approach of solving (30), (31) by iterations.We consider fork type symmetric bifurcate as shown in Fig. 4 where diameter of the parent branch is 2cm, whereas the peripheral branches are twice thinner.Elastic modulus is 1MPa, thickness 1.6mm.Boundary conditions have been specified as a flow at the inlet section to the system, (proximal flow, Fig. 4 ), and a pressure at each exit section, (distal pressure, Fig. 4).Calculations used 100 cell, Courant number CFL=2; ϴ=0.5.Calculated proximal pressure and a distal flow velocity, will be compared with in-vitro experimental testing data in future work.

Conclusion
Analytical and numerical solutions are presented and verified for the one-dimensional direct and inverse problems of fluid dynamics in an elastic thin-walled tube.Based on direct computational analysis the inverse solution was obtained, calculating the local elastic properties of the arterial wall using typical diagnostic measurements.
14) quantifies the correction to the classical Moens-Korteweg model, introduced by the present theory.While the linear model predicts all waves travel with the Moens-Korteweg speed ( MK c ), accounting for nonlinearities predicts speed of propagation which exceeds MK c by 10-15 % within the physiological range of transmural pressure.The speed -amplitude relationship for the peak of flow velocity u A , volume flow rate Q A , and pressure p A are retrieved from equations (

Figure 2 . 3 . 2 .
Figure 2. Numerical solution of Burgers equation.Comparison with an exact solution (a); Shock resolution (b)3.2.Flow in elastic arterial segmentsFig 3 shows dependence of a peak displacement amplitude on the flow velocity amplitudes based on the traveling wave model (15-17) and TVD approach.This solution serves as a test for the TVD model applied to the flow simulation inside elastic cylindrical cylinder.The downstream condition is specified as a non-reflected boundary condition, by setting to zero the backward moving Riemann invariant.On a left boundary the sinusoidal profile for the velocity was specified.We used Courant number CFL= 2, 100 cells along the tube, ϴ=0.5

Figure 3 .
Figure 3.A comparison between the Non-dimesionalized peak deflection magnitude and peak pressure magnitude the corresponding pressure and velocities values.The lower case letter variables are known at each step of iteration process.The total set of 6 equations to calculate 6 unknowns comprises three conservation laws for the flow rate and momentum ( L h -losses associated with the area jump and a flow direction change inside junction)

Figure 4 . 3 . 4 .
Figure 4. Boundary conditions specified as a flow at the inlet section (u1) to the system, and a pressure (p2,3) at each exit section3.4.Inverse problem of hemodynamicsUse of the model to solve the inverse problem of hemodynamic results in calculation of the local elastic modulus by non-invasively measuring pressure wave velocity, peak pressure, normal deflection or ejection time and ejection volume.To determine elastic modulus we use its contour plot in terms of PWV and ejection time built based on equations (12) or(13).For the specific values of an ejection time and a peak pressure we simply identify elastic properties with the nearest contour line, as shown in Fig.5a.All three combinations of PWV and ejection time result in the quantification for elastic modulus of about 1.04MPa as it is shown in Fig 5b.

Figure 5 .
Figure 5. Elastance surface plot from the model over a range of PWV and Ejection Time (a).Comparison of the model to three measured ejection times and PWV's (b)