The fluid mechanics involved in flow through a myocardial bridge is complex, because of the three dimensionality of the deformations, coupling of the fluid with the arterial wall and flow separation. To understand the complicated behaviour of the tube flow, it is convenient to start with a one-dimensional approximation, which qualitatively predicts the overall aspects in spatially averaged flow variables. It is commonly derived by using equations for mass and momentum conservation and can be found in several places [47–50]. Because the equations are in common use, we will shortly repeat the main assumptions made. However, we will point out the differences introduced by the wall deformation, the entrance type of flow and flow separation. Firstly, we assume that blood flow in reasonably large vessels can be modelled as incompressible, Newtonian fluid with constant density ρ
0 = 1055 kg/m3 and constant dynamic viscosity μ = 4 mPa s; the kinematic viscosity is defined as ν = μ/ρ
0. The Reynolds number Re
d
= ud/ν, based upon the vessel diameter, d and the mean flow velocity u, is below 2000 in non-diseased locations of the coronary arteries, so that the flow can be assumed to be laminar [9].
The basic geometry
Theories of longitudinal waves in tubes, with or without non-uniformities, non-linearity and frictional dissipation, are based on the idea that variation of the pressure in an axial cross-section is negligible. If the internal and external pressure along the main flow axis, x, of the artery at time t are given by p
int
(x, t) and p
ext
(x, t) respectively, the transmural pressure across the wall of a tube is p
tm
= p
int
- p
ext
. Henceforth we assume that the external pressure is constant in space and time and consequently it is the internal pressure whose gradients produce fluid acceleration. Our first geometrical simplification for modelling blood flow in arteries is that the curvature along the axis of the tube is assumed to be small everywhere and that the flow in the cardiovascular system is unidirectional, so that the problem can be defined in one space dimension along the x-axis. According to this we have simplified the anatomy of the myocardial bridge as shown in Figure 2.
The two arrows in Figure 2 denote the location of either circular (B - B) or oval (C - C) cross-section of the tube. Due to the fact that the wall thickness, h
0, is small compared to the bending radius R
d
(h
0
/R
d
≪ 1), we assume that the bending stress inside the wall is negligible. Consequently the cross-section of a circular tube (Figure 3 (left)) under deformation in z-direction, is given by the composition of a rectangle with two semicircular ends as illustrated in Figure 3 (right). This is consistent with the predominately eccentric deformation of bridged segments found in [3]. We note that negligence of bending stress causes the tube to collapse significantly earlier, i.e. the assumption is only satisfied if p
tm
≥ 0.
The deformation distance between the squeezing muscles and the breadth of the flat portion are denoted by D (x) and B (x, t) respectively. The equilibrium geometry of the cylindrical tube in Figure 3 (left) is characterised by the inner radius, R
0, the circumference U
0 = 2 π R
0 and the cross-section A
0 = π R0
2.
However the equilibrium cross-sectional area of the deformed tube is A
d
(x, t) (see Figure 3 (right)). The total cross-sectional area in the yz-plane of the tube is defined by A (x, t) = ∫
A
da and the actual circumference is U
p
(x, t) = 2(π R
d
+ B). Consequently the average flow velocity u (x, t) = 1/A ∫
A
ν
x
da, where ν
x
is the local value of the flow velocity in axial direction. The volume flux across a given section therefore is q (x, t) = A u.
As shown in the angiography 1 and Figure 2 the coronary arteries in myocardial bridges are structured by several wall deformations. Their number, degree and extension may independently vary with time, so that the axial curvature of the arterial wall for each of the n = 1...N myocardial bridges in series is characterised by N functions. The deformation is specified by a parameter ζ, defined as ζ = R
d
/R
0, which is chosen in the stenosis n to vary with time as
where g
n
(t) are periodic functions describing the temporal contraction of the muscle fibres. ζ
systole
and ζ
diastole
are fixed geometrical parameters between 0 and 1, specifying the degree of systolic and diastolic deformation respectively. We note that in the centre of the deformation ζ (x = x
s 2) = ζ
0 = R
d
/R
0, i.e. the degree of deformation increases with decreasing ζ
0 and consequently ζ
systole
<ζ
diastole
. To represent the time dependence of the deformation found in intra-vascular ultra-sound (IVUS) measurements [51], a synthetic deformation waveform g
n
(t) given by m = 1..3 sine/cosine harmonics is used.
Here Δt
n
is the time shift for each deformation with respect to the cardiac cycle (see Figure 4) and φ
m
are the phases in radian, chosen to be φ
1 = 3.5, φ
2 = 1.5, and φ
3 = 3.9. The axial curvature of each deformation is approximated by two hyperbolic tangent functions, so that
where x is the axial coordinate and x
tn
are the transition locations. Equation (3) smoothes the transition between the segmental domains Ω
n
by a transition length 1tn. The actual state of deformation can also be expressed by the ratio, ε, based upon the half-axes of the non-circular duct.
The pressure-area relationship
Due to its complicated structure, it is difficult to provide a synthetic mathematical description for the mechanical behaviour of vessel walls. Here, we focus on the most relevant structural features and the simplest mathematical model for arterial tissues. In the following we derive an algebraic pressure-area relationship for a vessel under external deformation. In this context the distensibility characterises the relative change in cross-sectional area with respect to the pressure increment for a given deformation according to A = A (ζ, p). If we assume that A' is the perturbation about the equilibrium area A
d
, the total cross-sectional area can be written as A (ζ, p) = A' (ζ, p) + A
d
(ζ). For a homogeneous, thin-walled (h
0
/R
0 ≪ 1), linear elastic tube, the stresses in circumferential direction are large compared to stresses in longitudinal direction and the hoop stress per unit length of the tube is
where E is the elastic modulus and σ is the Poisson ratio, which for practically incompressible biological tissue is approximately 0.5. The circumferential strain in the vessel wall is
Equation (5) can be rearranged into the form
whereas Equation (6) leads to an expression for the pressure dependence of the breadth, B, of the flat portion of the tube.
The total cross-sectional area is
A (R
d
, p) = π R2
d + 2 B (R
d
, p) R
d
. (9)
In the unperturbed state, the cross-section is
A
d
(R
d
) = R
d
U
0 - π R2
d. (10)
We can finally write the pressure induced perturbation as
It should be noted that under the assumption of linear elastic material with constant elastic modulus, equation (9) and (11) have the property that the area increases linearly with transmural pressure. Real arteries, however, resist over-expansion by having an incremental elastic modulus, E(ε), that increases with increasing strain [52]. It should be further noted that the area perturbation in equation (11) is not only dependent on pressure variation but also on the degree of deformation through R
d
. By using equation (7) and (9) we can finally write the pressure-area relation
Elastic properties of the coronary arteries
The elastic properties for a given section of the circular tube are obtained by using estimates for the volume compliance as suggested in [53], where the empirical approximation in exponential form is
In these estimates k
1, k
2, and k
3 are constants. With data for the volume compliance from Westerhof [54], Stergiopulos [55] and Segers [56] we obtain k
1 = 2.0 * 106 Kg s-2m-1, k
2 = -2.253 * 103 m-1, and k
3 = 8.65 * 104 Kg s-2m-1. A functional relationship for the wall thickness subject to the vessel radius can be found in [57], where
h
0 = a R0
b. (14)
The parameters for a = 3.87 and b = 0.63 were obtained by a logarithmic fit to data including vessel radii between 100 μm and 3000 μm. Equations (13) and (14) are used to determine the wall thickness and elastic properties of the vessel, if the radius is known. The assumption of small bending resistance is well satisfied if a R0
(b-1) ≪ 1, which for typical vessels under consideration is below 0.21.
Interaction of viscous boundary layer and inviscid core flow
In the following we investigate the solutions of the unsteady boundary layer equations by using an approximate integral method proposed by Veldman [58]. For this purpose, the potential flow of the two-dimensional equations governing the unsteady incompressible laminar boundary layer flow in axial symmetry is assumed to be in power-law form. By the introduction of similarity variables and the assumption that the evolution of the velocity profile is weakly dependent on x, the boundary layer equations reduce to the Falkner-Skan equation. Based on this ordinary differential equation, closed form solutions to the von Kármán integral momentum equation are obtained by a curve fit representation. The steady skin friction coefficients and the non-linear momentum correction coefficients corresponding to the velocity distributions are obtained and compared with known results. In particular the results of the steady solution are found to compare favourably with the Blasius solution and values for fully developed flow.
Boundary layer equations
The notion of the boundary-layer approximations was first developed by Ludwig Prandtl in the early 1900's. These well-known approximations [59] are applied widely in fluid mechanics. As the flow rate in the tube increases (i.e. high Reynolds number) the boundary-layer approximations become increasingly valid. The derivation of the axial boundary layer equations was first given by Mangler (1945) and can be found in [60]. In cylindrical coordinates (x, r, φ) with the corresponding velocity components (ν
x
, ν
r
, ν
φ
) and circumferential velocity ν
φ
= 0 (no swirl in S), they are
with ν
x
(x, R
d
, t) = 0 ν
r
(x, R
d
, t) = ν
w
(x, t), (18)
where R
d
(x, t) is the body shape along a xr-section of the tube or local surface radius measured from the axis and τ represents the shear stress, which is defined as τ = μ ∂VX/∂r. The derivation assumes that R
d
is much larger than the boundary layer thickness δ. The three-dimensionality of deformation generally makes it difficult to find a satisfactory solution for every compartment of the neither circular nor flat duct. However, considering severe deformation (e.g. ζ
0 = 0.2) the circumferential length of the flat portion of the vessel exceeds the circumferential length of the circular portion by a factor of four (
), thus we assume plane wedge flow for the calculation of viscous forces. Consequently ν
w
is taken to be the velocity component normal to the flat portion of the wall, which is ν
w
= -∂Rd/∂t. At the edge of the boundary layer, the free-stream velocity V (x, t) must be related to the pressure by the potential flow relation
Integral momentum equation
The integral momentum relation of von Kármán (1921) is obtained by multiplying the continuity equation (15) by u - V and subtracting it from the momentum equation (16). Integration over the bending radius and introduction of the integral relations for the displacement thickness
the momentum thickness
and the total displacement thickness δ** = δ* + θ
leads to
with u(0, t) = V(0, t), θ(0, t) = δ*(0, t) = 0, (24)
where c
f
is a non-dimensional friction factor defined as
The only difference to plane flow is the term involving ∂R
d
/∂
x
. If R
d
→ ∞ or ∂R
d
/∂x → 0, equation (23) reduces to the von Kármán integral momentum equation for plane flow. Compared to the frictional term cf/2 the influence of the term involving ∂R
d
/∂x on the boundary layer properties is indeed small (below 0.3%), thus disregarding the term for simplicity is appropriate. The boundary conditions in (24) assume a uniform inflow profile.
Falkner-Skan equation
Suitable solutions to the boundary layer equations in either plane or axial symmetry are found by the introduction of the stream function. A suitable coordinate transformation turns the equation for the stream function into the Görtler equation (derivation see [59]). The following approach mainly consists in assuming that the flow is locally self-similar and that it depends weekly on the coordinate x, so that the velocity profiles can be mapped onto each other by suitable scaling factors in y. Falkner and Skan have found a family of similarity solutions, where the free-stream velocity is of the power-law form
V(x) = Cxn, (26)
with a constant C and the power-law parameter n. The similarity variable η ~ y/δ(x) is set as
where f(η) is the dimensionless stream function and the prime refers to derivative with respect to η. The coordinate normal to the plate is denoted by y. However there are other general similarity solutions including the temporal dependence of the profile evolution [61]. The above similarity variables turn the boundary layer equations into a non-linear ordinary differential equation of order three, which is known as the Falkner-Skan-Equation
The parameter β is a measure of the pressure gradient ∂ p /∂ x. If β is positive, the pressure gradient is negative or favourable, β = 0 indicates no pressure gradient (i.e. the Blasius solution to flat plate flow) and negative β denotes a positive or unfavourable pressure gradient. We note that by assumption, β should vary slowly with coordinate x. The solutions are found numerically by a shooting method with f''(0) as free parameter by Hartree [62]. To avoid extensive calculations we follow a curve fit representation of three quantities extracted from solutions of the Falkner-Skan equation used in [19]:
The relation to flat plate flow is generally given by the shape factor, which is defined as H = δ*/0, whereas H
0 is the equivalent value for plane flow over a flat plate. The curve fits provide a good approximation for values of H between 1 and 20. At the separation point the wall shear stress vanishes, i.e. ∂VX/∂r = 0, which is equivalent to a shape factor H = 4. Relation (32) is not required for the calculations, however it is useful to predict the actual boundary layer thickness δ
99, where the fluid velocity differs by 1% from the free stream value. The relation for the shear stress given in [19] is
consequently the friction factor is
A schematic illustration of the actual flow profile along the tube axis is shown in Figure 5. We have chosen a uniform inflow profile with velocity V(0, t). The boundary layer (dashed line) grows from the leading edge, decreases in the converging part, while it grows in the divergent part of the tube. The upward triangles (▲) denote the point of separation, while downward triangles (▼) indicate the reattachment of the boundary layer. After separation the flow field can be seen as a top hat profile in the centre and a recirculation zone close to the walls. Due to the adjacent converging part the reattachment is forced early, because fluid is accelerated. In contrast the reattachment after the second diverging part takes place further downstream.
Averaged flow equations
The simultaneous viscid-inviscid boundary layer approach assumes an inviscid core flow, which follows equation (19) and a viscous boundary layer, which may be found by the solution of equation (23). The one-dimensional equations commonly used to simulate unsteady, incompressible blood flow in elastic tubes with frictional losses [53, 63] are given in averaged flow variables as
where F
ν
is the viscous friction term and χ is the momentum correction coefficient. The viscous friction term is defined as
and the momentum correction coefficient is
We rearrange the equations written in area and flow rate in terms of area and area-averaged axial flow velocity so that
The derivative of A
d
with respect to time in equation (39) is a prescribed function depending on R
d
(x, t). It is responsible for the volume displacement caused by the forced deformation of the tube.
Hagen-Poiseuille viscous friction and momentum correction
The determination of viscous friction factor and momentum correction coefficient requires knowledge about the velocity profile. For pulsatile laminar flow in small axially symmetric vessels a flow profile of the form
is used [50]. Here û is the free stream value of the axial velocity and R is the actual radius of tube, while γ is the profile exponent, which for a Hagen-Poiseuille flow profile is equal to two. Consequently the friction term is given by
F
ν
= -2 π ν(γ + 2) u = K
ν
u, (42)
whereas the friction coefficient for a parabolic profile is K
ν
= -8 π ν. The corresponding momentum correction coefficient is given by
which is 4/3 in the parabolic case. Such factors can be used for the purpose of correlating other variables as well as for direct calculation of pressure drop. We note that in the presence of a stenosis the total losses under the assumption of Hagen-Poiseuille flow are underestimated [55]. This comes mainly from underestimating the viscous forces and disregarding the losses caused by flow separation at the diverging end of the stenosis [64, 65].
Boundary layer derived viscous friction and momentum correction
Developing flow conditions in ducts of multiply connected cross-sections generally make it difficult to use the right friction factor. A variety of cross-sections are discussed in [45]. In those situations the similarity parameters are preferably based on the free stream velocity and the square root of the cross-sectional flow area as characteristic length scale i.e.
. Consequently we define the Reynolds number
, the Womersley number or frequency parameter
and the Strouhal number
. Here ω is the angular frequency defined as ω = 2πf, with f the base frequency of pulse wave oscillation. In other words we have multiplied Re and Sr by a factor of
, while Wo is multiplied by
. In the calculations we have given the Reynolds number inside the stenosis, Re
st
, based on
.
As previously mentioned the surface line of the flat portion of the non-circular duct dominates the circular portion at severe deformations, so that the computation of viscous forces is based on plane wedge flow. Consequently the thickness of the boundary layer is estimated in the xz-plane. Further we assume that the boundary layer has constant thickness along the circumference as illustrated in Figure 6. The latter assumption allows a simple derivation of the momentum correction coefficient and the viscous friction term. Integration over the cross-section leads to geometric relations for the areas occupied by the displacement thickness δ* and the total displacement thickness δ** in that cross-section. They are expressed as
A
δ*= 2B δ* + π [R2
d - (R
d
- δ*)2], (44)
A
δ**= 2B δ** + π [R2
d - (R
d
- δ**)2], (45)
It is obvious that A
d
> A
δ**and A
d
> A
δ*have to be satisfied to make sure that the flow is not fully developed. The momentum correction coefficient can be found by satisfying mass conservation for the mean flow and the core flow by
V (A - A
δ*
) = A u, (46)
which can be used together with equation (22) and (45) in the definition for the momentum correction (43), so that
The uniform inflow profile is identical to χ = 1, while the developing profile reaches its far downstream value of 1.39 after the entrance length within less than 4.5% from the analytical solution for the parabolic flow profile given in equation (43). We note that in the linearised system the total cross-section A in equation (47) is replaced by A
d
. According to equation (34) the friction factor built with the pressure dependent surface line, U
p
(x, t) is
Computations in a uniform tube show good agreement to the friction factor of the parabolic profile given in equation (42). After the entrance length the friction factor computed via the boundary layer theory reached its far downstream value to within 7%. Additionally the Fanning friction factor Reynolds number product in a deformed vessel geometry agrees well with experiments carried out in [45]. In contrast to other proposed models the underlying model does not require knowledge about the fully developed friction factor Reynolds product c
f
Re
fd
[45], or the incremental pressure drop factor K
∞ [66, 67]. The above results for momentum correction and viscous friction quantify the entrance conditions typically encountered in studies of the arterial system.
Validity
The stationary boundary layer approximation becomes increasingly valid if the Reynolds number increases and when the ratio of unsteady forces to viscous forces given by the Womersley number is small. In the left coronary artery the values for Womersley and Strouhal number, built with pulsatile frequency were around 4 and below 0.2 respectively. The approximation of the actual flow profiles by near equilibrium flow profiles is justified for stationary flow, however, for time dependent flow the period T of the deformation function in equation (3) has with be large compared to the viscous diffusion time through the boundary layer: t
d
= δ2/V ≪ T. This is well satisfied for the situation under consideration, where t
d
in the centre of the deformation was between 4.4 * 10-3 s and 0.2 s at values of ζ
0 of 0.25 and 1 respectively.
Computational domain
Based on 83 angiographies, Dodge et al. [68, 69] presented a normal anatomic distribution of coronary artery segments and proposed a terminology, which we used for our model of the left coronary artery (LCA): the left main coronary artery (LMCA) bifurcates into the left anterior descending artery (LAD) and the left circumflex artery (LCxA). The main branches of the LAD include the 1st, 2ndand 3rddiagonal branch (Dl, D2, D3) and the 1st, 2ndand 3rdseptal branch (S1, S2, S3). The main branches of the LCxA include the 1stand 2ndobtuse marginal branches (OM1, OM2). The exact intrathoratic location and course of each one of the 27 arterial segments and branches of the LCA are illustrated in Figure 7. We note that in the present 1D approximation the arterial tree is composed of tubular entities. The branching angles in figure 7 serve only artistic purposes.
Boundary conditions
Because the coronary flow is primarily driven by the aortic pressure, the pulsatile inflow condition to the LMCA was represented by a periodic extension to a synthetic pressure wave in the exponential form
where p
s
is the static pressure and p
0 is the amplitude of the exponential waveform, while t
r
is the rising time. We have chosen the parameters according to measurements in [38]. The baseline condition is represented by a stationary pressure p
s
= 8 kPa and a pressure amplitude of p
0 = 7 kPa, while for the inlet pressure under dobutamine challenge we have chosen p
s
= 6.6 kPa and p
0 = 5.5 kPa. In both cases the raising time was t
r
= 0.25 s. The pressure wave at the inlet is shown in Figure 8.
The branching conditions between the 1D entities are implemented by the requirement of constant pressure at the branching point and mass conservation throughout the bifurcation [63]. To avoid wave reflections at the ends of the tubes the boundary conditions are implemented by a characteristic system of one-way wave equations [46]. There are several ways to account for peripheral reflections at the terminals of a vascular network. In the current simulations we have implemented a three-element windkessel model for the termination of the left coronary arterial tree [70, 71]. The main advantage of this model is to consider the compliant-capacitive effects due to micro-vessels and arterioles. To satisfy the Blasius solution at the leading edge of the tube, we have assumed a uniform flow profile at the entrance (V(0, t) = u(0, t)).
Numerical implementation
Due to the non-linear terms in equation (23) and (40) the solutions for haemodynamically developing flows are generally more difficult to obtain than fully developed flows or oscillating flows with a frequency dependent Stokes boundary layer. Developing flows require simultaneous solution of the momentum equation (39), the continuity equation (40) and the integral momentum equation (23), together with the boundary conditions given in (49) and (24) respectively. The system of equations cannot be solved analytically, so that the interior domain was solved by a second order predictor-corrector MacCormack finite difference scheme with alternating direction for prediction and correction in each time step [72]. To implement the boundary and interface conditions it is convenient to disregard viscous friction and rewrite equation (39) and (40) in terms of characteristic variables [46]. The momentum correction factor in equation (47) and the viscous friction in equation (48) are given by the solution to the integral momentum equation (23) and the two curve fits to the Falkner-Skan equation in (30) and (31). They are solved by discretisation using the same second order MacCormack scheme and iterative solution of the resulting set of discrete non-linear equations by a combined root bracketing, interval bisection and inverse quadratic interpolation method of van Wijngaarden-Brent-Dekker. To start the computation the Blasius solution at each time step provides values for the boundary layer thickness a few grid points downstream of the entrance. The solution was applied to the steady integral momentum equation as boundary condition for δ*. Downstream marching the solution leads to the values of boundary properties along the tube axis, which through equation (47) and (48) provide the required values of momentum correction and viscous friction to the averaged flow equations respectively.