Unsteady Integrodifferential Equation of Fluid-Structure Interaction in Constricted Collapsible Tube Model of Diseased Human Coronary Artery

Abstract

Unsteady flow in a collapsible tube is analyzed to simulate a diseased human coronary artery. The novelty of the approach is that the set of equations governing the fluid-structure interaction is reduced to a single integrodifferential equation in the transient state. The equation is then solved using the finite difference method to obtain the flow characteristics and compliant wall behavior. Three control parameters are investigated, namely, Reynolds number, inlet transmural pressure, and the wall thickness. The predicted wall deflection is quite large at low Reynolds numbers, suggesting possible approach to breakdown in equilibrium. The transmural pressure increases with wall deflection and bulges appear at the ends of the membrane indicating critical stage of stability, consistent with previous studies. Increase in wall thickness reduces the wall deflection and ultimately results in its collapse which may indicate another breakdown in equilibrium. An increase in internal pressure is required to maintain membrane stability.

1. Introduction

This paper describes a two-dimensional (2D) analytical study of collapsible tube simulating unsteady flow in a stenosed coronary artery. The effect is investigated of peristaltic flow in a geometrically nonlinear elastic tube, whose walls deform due to the transient transmural pressure. The arterial wall motion results in a strong fluid-structure interaction. The study is based on the hypothesis that the transient fluid-structure interaction in such a channel could be adequately represented as a single integrodifferential equation which can be solved using finite-difference numerical method. The novelty of the study is the analysis of unsteady flow and the methodology employed to solve the resulting integrodifferential equation.

The problem posed by fluid flow through flexible tubes with relatively thin walls is theoretically challenging and practically significant. The theoretical challenges are due to the interaction between the fluid flow and the elastic channel walls. Thus, fluid flow and structural parameters have to be computed simultaneously. In addition, the boundary conditions cannot be completely defined in advance, due to the continuously evolving boundary. The practical significance of these problems are considerable, in particular in bioengineering and biomedical systems because of their pivotal role in describing fluid flow in living organs including cardiovascular, respiratory, and urinary systems [1]. In particular, vessels experiencing compressive transmural pressure and consequent remodeling include the vein above the heart, the coronary arteries during systole, and the airways during forced expiration [2]. The early theoretical investigations [3, 4] despite reproducing flow characteristics were unable to predict the details of wall deformation. Most previous studies have also focused on steady flow situations in order to reduce the computational complexity. The present paper relaxes these simplifications by analyzing the unsteady flow of deformable tube and solving the resulting integrodifferntial equations.

Flows in collapsible tubes have been thoroughly reviewed in a previous study [5]. The approach used in solving such flows depends on the flow situation. Steady flow with low Reynolds number and negligible inertia is referred to as lubrication-type flow [6]. Unsteady flow in which variations in the flow direction are significant is called Stokes-type flow [6]. In this case, the Reynolds number is typically high and inertia is taken into consideration. The full Navier-Stokes equations must be solved in such a situation, requiring numerical analysis. The elastic channel wall is usually modeled as a membrane using the thin-shell theory with negligible bending stiffness.

Several methods of modeling flow in flexible tubes have been summarized in a recent study [1]. The study discusses the significance of the interaction between fluid flow within the tube and the wall deformation. The physiological implications of such fluid-structure interaction are also fully explored, including human cardiovascular and respiratory systems. The study concludes that despite the intense interest, the diverse mechanisms generating instabilities in single-phase flow through a flexible tube remain only partially understood. Therefore, there is a need for more systematic investigations that will shed light on generic relationships and experimental observations as well as physiological applications carefully accounting for the mechanical properties of tissues.

Modeling of flow in collapsible tubes has ranged from one-dimensional (1D) to 3D [7–17]. Shapiro [7] investigated 1D steady flow in a partially collapsed thin-walled tube and demonstrated that steady flow cannot exist in a collapsible tube under flow velocity as high as the speed of propagation of small amplitude pressure waves. Under these conditions, choking and unsteadiness set in. Two-dimensional studies [8–11] include zero Reynolds number steady flow or finite Reynolds number unsteady flow, in which part of the rigid wall was replaced by an elastic membrane. Pedley [12] investigated the breakdown steady flow using lubrication theory in a symmetric 2D channel. The resulting nonlinear ordinary differential equations were solved numerically. The results suggest that reducing the longitudinal tension in the channel wall increases the deformation for the same flow rate and external pressure. Djordjevic and Vukobratovic [13] studied of steady, viscous flow in a collapsible channel at both low and high Reynolds numbers. The elastic channel wall behavior was modeled by a geometrically nonlinear Karman shell theory. In each case, a single integrodifferential equation was derived to describe the wall configuration and the governing equation was numerically solved. The wall configuration was predicted as a function of some control parameters. Other studies considered steady deformations of the fully three-dimensional system [14, 15]. The tube wall was modeled as a circular cylindrical shell and geometrically nonlinear shell theory was used to model its large nonaxisymmetric deformation (post buckling). The fluid flow was modeled using lubrication theory assuming low Reynolds numbers and a small wall slope in the direction of the flow. Pedley [16] investigated flow through collapsible arteries. More recently, perturbation and numerical methods were used to examine the drag reduction in gel-lined tubes and other surfaces [17].

Most of the variety of models applied to collapsible tube have used as the pivotal parameter, the transmural pressure, defined as the difference between the external and internal pressure. If this pressure differential is negative the channel experiences internal pressure, giving it a spherical cross-sectional shape [18]. If the pressure differential is positive the channel walls collapse. Many of the previous studies [7, 12, 16] have also focused on steady flow problems due to the complexity of the equations governing the fluid flow and wall configuration, and in order to reduce the number of control parameters. Other studies considering unsteady flow problems employed the Finite Element Method (FEM) to solve the governing equations [1, 10]. Application of FEM to the type of problem considered is complicated by the variety of competing factors including interpolations at the nodes of discretized elements, strain-displacement relation, stress-strain calculation, element equilibrium equations, global equilibrium equations, and so forth [19].

In the present study, unsteady 2D Stokes equations will be derived and used to analyze the fluid-structure interaction at relatively low Reynolds number. The governing equations reduce to a single integrodifferential equation which is solved numerically by the finite-difference method (FDM). The advantage of FDM is its recursive approach which allows convenient and accurate representation of the nonlinear fluid flow equations. The effect on flow and wall configurations are investigated of three control parameters, the Reynolds number (or equivalently, the inlet velocity), transmural pressure, and wall thickness.

2. Formulation

2.1. Hypothesis

The transient fluid-structure interaction in a collapsible tube model of diseased human coronary artery could be represented as a single integrodifferential equation which can be solved using finite-difference numerical method.

This section describes the derivation and solution of the governing equation. The transmural pressure distribution on the channel wall is determined from lubrication theory. The nonlinear Karman plate theory is used to express the wall elastic behavior. A single integrodifferential equation is derived for the first time, governing the membrane behavior in the transient state. The equations for pressure and compliant wall are solved using finite difference methods. This analytical approach differs from the relatively expensive application of finite-element method as in previous studies.

2.2. Geometry

Consider an idealized symmetrical diseased human coronary artery (Figures 1 and 2) through which an incompressible, viscous fluid is flowing. This two-dimensional (2D) model is inspired from Djordjevic and Vukobratovic’s [13] where the flow is assumed to be fully developed at the entrance (𝑥^∗=0) of the segment and remains laminar throughout. Thus the velocity profile in the segment is independent of Reynolds number. The artery wall configuration is driven by stresses acting on the surface. The elastic wall configuration is assumed to be a plate. Thus both wall inertia and fluid inertia are represented in the calculation.

Figure 1 

Schematic sketch of 2D collapsible channel.

Figure 2 

Segment of channel wall.

We employ lubrication theory, assuming small wall slope and negligible wall inertia [12, 13]. The segment has length 𝑙^∗ comprising parallel compliant walls with elasticity modulus 𝐸. In Figure 1, 𝑥^∗and𝑦^∗ are the Cartesian coordinates, 𝑢^∗(𝑥^∗,𝑦^∗,𝑡^∗)and 𝑣^∗(𝑥^∗,𝑦^∗,𝑡^∗)are the velocity components in 𝑥^∗and𝑦^∗ directions, respectively, and 𝑝^∗(𝑥^∗,𝑦^∗,𝑡^∗) is the pressure. The parameters𝜌, 𝜈, and 𝑝_𝑎are the density, kinematic viscosity, and ambient or external pressure, respectively. The parameter 𝑝_𝑜 is the constant entrance pressure at 𝑥^∗=0, 𝑢_𝑜 the average entrance velocity, 𝛿^∗(𝑥^∗,𝑡^∗) the arbitrary position of the channel wall with respect to direction 𝑥^∗and time 𝑡^∗, and 𝛼(𝑥^∗,𝑡^∗) is the angle between the bending wall and the principal axis 𝑥^∗. The asterisks (*) denote dimensional parameters that will be nondimensionalized later.

2.3. Governing Equations: Fluid Flow

The incompressible pulsatile flow is modeled using time-dependent, Navier-Stokes equations for Newtonian fluids coupled with the continuity equation. The equations are nondimensionalized using 𝑢_𝑜,𝛿_𝑜,𝑝_𝑜, and 𝑡_𝑜as scales for velocities, lengths, pressure, and time, respectively. The nondimensional equation can be expressed thus: 1√𝜀𝜕𝑢𝜕𝑡+𝑢𝜕𝑢𝜕𝑥+𝑣𝜕𝑢𝜕𝑦=−Eu𝜕𝑝+1𝜕𝑥𝜕Re²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²,1√𝜀𝜕𝑣𝜕𝑡+𝑢𝜕𝑣𝜕𝑥+𝑣𝜕𝑣𝜕𝑦=−Eu𝜕𝑝+1𝜕𝑦𝜕Re²𝑣𝜕𝑥²+𝜕²𝑣𝜕𝑦²,𝜕𝑢+𝜕𝑥𝜕𝑣𝜕𝑦=0, (2.1) where the nondimensional parameters 𝑥=𝑥^∗/𝛿_𝑜, 𝑦=𝑦^∗/𝛿_𝑜, 𝑢=𝑢^∗/𝑢_𝑜, 𝑣=𝑣^∗/𝑢_𝑜, 𝑡=𝑡^∗/𝑡_𝑜√𝜀, 𝑡_𝑜=𝛿_𝑜/𝑢_𝑜is the time scale, Eu=𝑝_𝑜/𝜌𝑢_𝑜², and Re=𝑢_𝑜𝛿_𝑜/𝜈 are Euler number and Reynolds number, respectively, and 𝛿_𝑜/𝑙^∗≪1, 𝛼_max≪1. The transmural pressure 𝑝 is 𝑝=𝑝𝑥,𝑦,𝑡^𝑎−𝑝^∗𝑥^∗,𝑦^∗,𝑡^∗𝑝_𝑜, (2.2) where 𝑝_𝑜 and 𝑝_𝑎 are the reference and atmospheric pressures, respectively.

The above Navier-Stokes equations will be simplified using the following assumptions: the stenosed segment of the channel is relatively long, (𝛿_𝑜/𝑙^∗≪1), and the maximum angle of the slope of the stenosis with respect to the 𝑥^∗-axis 𝛼_max≪1. The latter two conditions imply gradual change of the internal diameter of the channel along the principal axis. The relatively small inclination of the stenosis implies that all the physical quantities related to the structure will also be subject to very small variation in the direction of the principal axis. The flow is assumed to be fully developed throughout the channel. Thus the velocity component in the flow direction is large compared to the traverse component.

We translate small variations by introducing a coordinate 𝜉=𝜀𝑥and write the transverse velocity component in the form𝑣=𝜀𝑉(𝜉,𝑦,𝑡), where 𝑉=𝑂(1) and 0<𝜀≪1 will be defined later. We also assume that 𝛿_𝑜/𝑙^∗=𝑂(𝜀) and𝛼_max=𝑂(𝜀),𝑝(𝑥,𝑦,𝑡), the space and time-varying transmural pressure (outer atmospheric pressure minus inner fluid pressure) is much smaller than the reference pressure 𝑝_𝑜 which is on the order10⁵Pa [20]. We assume that the transmural pressure is transmitted through the wall alone. This condition can be formulated as 𝑝=𝜀^𝑛𝑃(𝜉,𝑦,𝑡),𝑃=𝑂(1), where the parameter 𝑛>1 will be defined later.

In terms of the new parameters, the Navier-stokes and continuity equations become 1√𝜀𝜕𝑢𝑢𝜕𝑡+𝜀𝜕𝑢𝜕𝜉+𝑉𝜕𝑢𝜕𝑦=−𝜀^𝑛+1Eu𝜕𝑃+1𝜕𝜉𝜀Re²𝜕²𝑢𝜕𝜉²+𝜕²𝑢𝜕𝑦²,𝜀√𝜀𝜕𝑉𝜕𝑡+𝜀²𝑢𝜕𝑉𝜕𝜉+𝑉𝜕𝑉𝜕𝑦=−𝜀^𝑛Eu𝜕𝑃+𝜀𝜕𝑦𝜀Re²𝜕²𝑉𝜕𝜉²+𝜕²𝑉𝜕𝑦²,𝜕𝑢+𝜕𝜉𝜕𝑉𝜕𝑦=0. (2.3) We consider the flow to be lubrication type, that is, the fluid is flowing very slowly and frictional forces are larger than inertial forces. The latter assumption implies that Reynolds number is low. In order to exploit the lubrication approximation, we make the following substitutions into (2.3): 1Eu=𝜀^𝑛+1,1𝜀Re=𝜆=𝑂^𝑘, (2.4) where 0<𝑘<1 and Re=𝑂(1).

Keeping terms of up to order 𝜀², we obtain 1√𝜀𝜕𝑢𝑢𝜕𝑡+𝜀𝜕𝑢𝜕𝜉+𝑉𝜕𝑢𝜕𝑦=−𝜕𝑃𝜕𝜕𝜉+𝜆²𝑢𝜕𝑦²+𝑂𝜀²,√𝜀𝜕𝑉1𝜕𝑡=−𝜀𝜕𝑃𝜀𝜕𝑦+𝑂²,𝜕𝑢+𝜕𝜉𝜕𝑉𝜕𝑦=0. (2.5)(2.6)(2.7) Equations (2.5), (2.6), and (2.7) have to be solved subject to the boundary conditions 𝑦=𝛿(𝑡,𝜉)∶𝑢=0,𝑉=𝜕𝛿^∗𝜕𝑡^∗=1𝜀^3/2𝜕𝛿,𝜕𝑡𝑦=0∶𝜕𝑢=𝜕𝜕𝑦³𝑢𝜕𝑦³=𝜕⁴𝑢𝜕𝑦⁴=⋯=0,𝑉=0. (2.8) Our main goal is to obtain the nondimensional pressure distribution 𝑃 in the channel wall. In order to achieve this goal, we first solve (2.5) for 𝑃 thus𝑃𝑡,𝜉,𝑦=𝑃𝑡,𝑦−^𝜉₀𝜆𝜕²𝑢𝜕𝑦²−1√𝜀𝜕𝑢𝑢𝜕𝑡−𝜀𝜕𝑢𝜕𝜉+𝑉𝜕𝑢𝜕𝑦𝑑𝜍, (2.9) where 𝑃(𝑡,𝑦)is an arbitrary function of integration.

Taking into account the small ratio of radius to length and the parabolic profile along the flow direction, we assume that fluid velocity in the horizontal direction is given by 𝑢=3𝜉,𝑦,𝑡2𝑈𝑦1−²𝛿², (2.10) where 𝑈(𝑡,𝜉) is the cross-sectional average velocity. The latter is calculated as follows. First, we integrate the continuity equation (2.7) from 0 to 𝑦and obtain 𝑉3𝑡,𝜉,𝑦=−2𝜕𝑈𝑦𝜕𝜉𝑦−³3𝛿²−𝑦³𝛿³𝜕𝛿𝜕𝜉𝑈+𝐶(𝑡,𝜉), (2.11) where 𝐶(𝑡,𝜉) is an arbitrary function of integration.

By imposing the relevant boundary conditions on 𝑉 that is (2.8), we obtain 𝐶()1𝑡,𝜉=0,𝜀^3/2𝜕𝛿+𝜕𝑡𝜕𝑈𝜕𝜉𝛿+𝜕𝛿𝜕𝜉𝑈=0. (2.12) Thus 1𝑈=𝛿𝜃(𝑡)−1𝜀^3/2𝜉₀𝛿_𝑡(𝑡,𝑧)𝑑𝑧, (2.13) where θ is an arbitrary function of integration and 𝛿_𝑡 denotes 𝜕𝛿/𝜕𝑡.

Next we substitute (2.10), (2.11), and (2.13) back into (2.9) giving 𝑃𝑡,𝜉,𝑦=𝑃𝑡,𝑦−^𝜉₀𝜆𝑦−3𝑈𝛿²−1√𝜀32𝜕𝑈𝑦𝜕𝑡1−²𝛿²𝑦+3²𝛿³𝜕𝛿𝑈3𝜕𝑡−𝜀2𝑈𝑦1−²𝛿²32𝜕𝑈𝑦𝜕𝑡1−²𝛿²𝑦+3²𝛿³𝜕𝛿𝑈+−3𝜕𝑡2𝜕𝑈𝑦𝜕𝜉𝑦−³3𝛿²−𝑦³𝛿³𝜕𝛿𝑈𝑦𝜕𝜉−3𝑈𝛿²𝑑𝜍. (2.14) In order to determine𝑃(𝑡,𝑦), we solve (2.6) for 𝑃. Then, we substitute the 𝑃 obtained into (2.14) and obtain after some manipulations, 𝑃𝑡,𝑦=^𝜉₀𝜆𝑦−3𝑈𝛿²−1√𝜀32𝜕𝑈𝑦𝜕𝑡1−²𝛿²𝑦+3²𝛿³𝜕𝛿𝑈3𝜕𝑡−𝜀2𝑈𝑦1−²𝛿²32𝜕𝑈𝑦𝜕𝑡1−²𝛿²𝑦+3²𝛿³𝜕𝛿𝑈+−3𝜕𝑡2𝜕𝑈𝑦𝜕𝜉𝑦−³3𝛿²−𝑦³𝛿³𝜕𝛿𝑈𝑦𝜕𝜉−3𝑈𝛿²𝑑𝜍−𝜀^3/2𝑦₀𝜕𝑉𝜕𝑡𝑑𝑦^′+𝑃̆(𝑡), (2.15) where 𝑃̆(𝑡) is an arbitrary function of 𝑡.

Finally, by substituting (2.15) back into (2.14) and solving the resulting equation, we obtain the general form of the pressure distribution written as 𝑃=𝑃̆(𝑡)𝑡,𝜉,𝑦+𝜀^3/232𝜕²𝑈𝑦𝜕𝑡𝜕𝜉2²−𝑦⁴12𝛿²+𝑦⁴4𝛿³𝜕𝛿𝜕𝑡𝜕𝑈−3𝜕𝜉4𝑦⁴𝛿⁴𝜕𝛿𝜕𝑡𝜕𝛿𝑦𝜕𝜉𝑈+⁴4𝛿³𝜕²𝛿𝑦𝜕𝑡𝜕𝜉𝑈+⁴4𝛿³𝜕𝑈𝜕𝑡𝜕𝛿.𝜕𝜉 (2.16) The nondimensional pressure distribution over the stenosed surface (𝑦=𝛿) is then 𝑃=𝑃̆(𝑡)𝑡,𝜉,𝑦=𝛿+𝜀^3/258𝜕²𝑈𝛿𝜕𝑡𝜕𝜉²+𝛿4𝜕𝛿𝜕𝑡𝜕𝑈−3𝜕𝜉4𝜕𝛿𝜕𝑡𝜕𝛿𝛿𝜕𝜉𝑈+4𝜕²𝛿𝛿𝜕𝑡𝜕𝜉𝑈+4𝜕𝑈𝜕𝑡𝜕𝛿,𝜕𝜉 (2.17) where 𝑈is given by (2.13) and 𝑃̆(𝑡) is an arbitrary function.

By substituting (2.13) into (2.17) we finally obtain the pressure distribution over the stenosed surface as 𝑃()=𝑃̆(𝑡)+1𝑡,𝜉,𝛿4𝜀^3/252𝜕²1𝜕𝑡𝜕𝜉𝛿𝜃(𝑡)−1𝜀^3/2𝜉₀𝛿_𝑡(𝑡,𝑧)𝑑𝑧𝛿²+𝛿𝜕𝛿𝜕𝜕𝑡1𝜕𝜉𝛿𝜃(𝑡)−1𝜀^3/2𝜉₀𝛿_𝑡()−3𝑡,𝑧𝑑𝑧𝛿𝜕𝛿𝜕𝜉𝜕𝛿𝜃(𝑡)−1𝜕𝑡𝜀^3/2𝜉₀𝛿_𝑡()+𝜕𝑡,𝑧𝑑𝑧₂𝛿𝜃(𝑡)−1𝜕𝑡𝜕𝜉𝜀^3/2𝜉₀𝛿_𝑡(𝑡,𝑧)𝑑𝑧+𝛿𝜕𝛿𝜕𝜕𝜉1𝜕𝑡𝛿𝜃(𝑡)−1𝜀^3/2𝜉₀𝛿_𝑡(𝑡,𝑧)𝑑𝑧. (2.18) Note that pressure is evaluated at points which happen to be also the limits of integration. Therefore, care must be taken when choosing the numerical method to evaluate (2.18). Here, we define a new variable (𝑡)−1𝐵=𝜃𝜀^3/2𝜉₀𝛿_𝑡(𝑡,𝑧)𝑑𝑧 (2.19) such that 𝐵 is completely defined by 𝜕𝐵(𝑡,𝜉)1𝜕𝜉=−𝜀^3/2𝜕𝛿(𝑡,𝑧),𝐵()(𝑡).𝜕𝑡𝑡,0=𝜃 (2.20)

2.4. Governing Equations: Artery Wall

The artery wall configuration is driven by hemodynamic stresses acting on its surface. The elastic wall configuration is solved as a plate. We start by applying Newton’s second law of motion to the surface element presented in Figure 2. It should be remarked that α in this section represents the angle of wall inclination as illustrated in Figure 2. Newton’s second law is expressed by^⇀𝐹=𝑀⋅𝑉^⇀̇, (2.21) where ^⇀𝐹 is the external force vector, and 𝑀 and 𝑉̇^⇀are mass and acceleration vector, respectively.

The transverse component of (2.21) is obtained through force balance on artery wall material element (see Figure 1) in the transverse direction, thus 𝜏_𝑤sin(𝛼)𝑑𝑥^∗−(𝐹+𝑑𝐹)sin(𝛼+𝑑𝛼)+𝐹sin(𝛼)−𝑝_𝑎−𝑝^∗𝑑𝑥^∗+𝑉+𝑑𝑉−𝑉=ℎ𝜌𝑑𝑥^∗𝑑²𝑏^∗𝑑𝑡^∗2, (2.22) where 𝜏_𝑤 is the shear stress over the inner face of the wall, 𝐹 is the magnitude of the force at inclination 𝛼,𝐹+𝑑𝐹 is the magnitude of the force at inclination 𝛼+𝑑𝛼 as a result of the bend, and 𝜌 the material density 𝐹=𝜎_𝑥∗ℎ,𝑉=−𝑑𝑀𝑑𝑥^∗,𝑀=^ℎ/2_−ℎ/2𝜎_𝑥∗𝑦^∗𝑑𝑦^∗,𝜎_𝑥∗=𝐸1−𝜈²𝜀_𝑥∗. (2.23) Based on Karman’s nonlinear shell theory [21], Djordjevic and Vukobratovic [13] have shown that the finite strain component in the axial direction takes the form: 𝜀_𝑥′=𝑑𝑎^∗𝑑𝑥^∗+12𝑑𝑎^∗𝑑𝑥^∗2+𝑑𝑏^∗𝑑𝑥^∗2. (2.24) We assume that 𝛼 is small. Hencetan𝛼=−𝑑𝑏^∗/𝑑𝑥^∗,sin(𝛼)≅𝛼 andsin(𝛼+𝑑𝛼)≅𝛼+𝑑𝛼. The latter approximation converts (2.22) into −𝜏_𝑤𝑑𝑏^∗𝑑𝑥^∗+𝐸ℎ1−𝜈²𝑑²𝑏^∗𝑑𝑥^∗2𝜀_𝑥∗−𝐸ℎ³121−𝜈²𝑑⁴𝑏^∗𝑑𝑥^∗4−𝑝_𝑎−𝑝^∗𝑑=ℎ𝜌²𝑏^∗𝑑𝑡^∗2+h.o.t, (2.25) where the abbreviation h.o.t stands for Higher Order Terms. We may recast (2.25) into 𝑇𝑑²𝑏^∗𝑑𝑥^∗2𝑑−𝐷⁴𝑏^∗𝑑𝑥^∗4−𝜏_𝑤𝑑𝑏^∗𝑑𝑥^∗−𝑝_𝑎−𝑝^∗𝑑=ℎ𝜌²𝑏^∗𝑑𝑡^∗2+h.o.t, (2.26) where 𝐷=𝐸ℎ³/12(1−𝜈²) and𝑇=(𝐸ℎ/(1−𝜈²))𝜀_𝑥∗ are the bending stiffness and the longitudinal tension, respectively.

The axial component (𝑋) of the equation of motion is similarly obtained through momentum balance on artery wall element in the axial direction, −𝐹cos(𝛼)+(𝐹+𝑑𝐹)cos(𝛼+𝑑𝛼)−𝜏_𝑤cos(𝛼)𝑑𝑥^∗+𝑝_𝑎−𝑝^∗(𝛼)𝑑sin=ℎ𝜌²𝑎^∗𝑑𝑡^∗2𝑑𝑥^∗+h.o.t. (2.27) Using the fact that 𝛼 is small we obtain 𝑑𝑇𝑑𝑥^∗−𝜏_𝑤+𝑝_𝑎−𝑝^∗𝜕𝑏^∗𝜕𝑥^∗𝑑=𝜌ℎ²𝑎^∗𝑑𝑡^∗2. (2.28) Now, we proceed to nondimensionalize equations (2.25) and (2.26) by introducing the following nondimensional variables: 𝑥𝑥=^∗𝛿_𝑜𝑎,𝑎=^∗𝛿_𝑜𝑏,𝑏=^∗𝛿_𝑜()1,𝑎=𝜀𝐴𝑡,𝜉,𝜉=𝜀𝑥,𝑡=√𝜀𝑡^∗𝑡_𝑜,𝑡_𝑜=𝛿_𝑜𝑢_𝑜. (2.29) In terms of the new variables, we have𝜀_𝜉=𝜀²𝑑𝐴+1𝑑𝜉2𝜕𝑏𝜕𝜉²+𝑂𝜀⁴. (2.30) Equations (2.25) and (2.26) become, respectively, 𝜀⁴𝐸1−𝜗²ℎ𝛿_𝑜𝑑𝐴+1𝑑𝜉2𝜕𝑏𝜕𝜉²𝜕²𝑏𝜕𝜉²−𝐸𝜀⁴121−𝜗²ℎ𝛿_𝑜³𝜕⁴𝑏𝜕𝜉⁴−𝜀^𝑛𝑃𝑡,𝜉,𝑦𝑝_𝑜=𝜀^𝑛𝑝_𝑜ℎ𝛿_𝑜𝑑²𝑏𝑑𝑡²+𝑂𝜀^𝑛+1, (2.31)𝜀³𝐸1−𝜗²ℎ𝛿_𝑜𝜕𝜕𝜉𝑑𝐴+1𝑑𝜉2𝜕𝑏𝜕𝜉²+𝑂2𝑝_𝑜𝜀^𝑛+1ℎ𝛿_𝑜=0. (2.32) Integration of (2.32) yields 𝑑𝐴+1𝑑𝜉2𝜕𝑏𝜕𝜉²=𝐾(𝑡). (2.33) Thus ^𝐿₀𝜕𝑏𝜕𝜉²𝑑𝜉=2𝐾(𝑡)𝐿, (2.34) where 𝐿 is the nondimensional wall length.

Using the dominant balance in (2.31), we set 𝜀⁴𝐸1−𝜗²ℎ𝛿_𝑜=𝜀^𝑛𝑝_𝑜. (2.35) Recall also that 𝑝Eu=⁰𝛿²₀𝜆²𝜌𝜈²=𝑂(1). (2.36) Equations (2.35) and (2.36) imply that 𝑛=−1+lnEuln𝜀,𝜀=𝜌𝜈²(1−𝜗²)𝐸ℎ𝛿₀𝜆^21/5. (2.37) Using (2.34), (2.22) can be rewritten as 𝐾𝜕²𝑏𝜕𝜉²−1ℎ12𝛿_𝑜²𝜕⁴𝑏𝜕𝜉⁴ℎ−𝑃=𝛿_𝑜𝑑²𝑏𝑑𝑡². (2.38) Recalling that 𝑏=𝛿−1, (2.38) may be transformed further to 𝐾𝜕²𝛿𝜕𝜉²−1ℎ12𝛿_𝑜²𝜕⁴𝛿𝜕𝜉⁴ℎ−𝑃=𝛿_𝑜𝜕²𝛿𝜕𝑡², (2.39) or equivalently −𝑚₁²𝛿^(IV)+𝑚₂𝛿^(II)−𝑚₃=12𝑚₁𝛿̈, (2.40) where the superscripts and the overdot stand for differentiation with respect to 𝜉 and 𝑡, and 𝑚₁=ℎ𝛿_𝑜,𝑚₂=12𝐾,𝑚₃=12𝑃. (2.41) In the next section, we describe the numerical methods used to solve the above governing equations.

2.5. Numerical Solution Procedure

The main issue with the wall governing equation (2.41) in the transient state is solving a single differential equation with two unknowns, including transmural pressure distribution 𝑃 and wall deflection𝛿. Computation of the wall shape (wall deflection) at the next time step proceeds as follows. Knowing the wall shape at the initial time step (deflection distribution 𝛿at the initial time), we solve the equation of the transmural pressure equation (2.18) for 𝑃. Then, at the next time step we solve the wall governing equation (2.39) for 𝛿 under the assumption that the time step Δ𝑡 is small enough to apply 𝑃|_𝑡=𝑛≈𝑃|_𝑡=𝑛+1. Next, given the updated solution for 𝛿, we can find the corresponding updated 𝑃, and so on until the final iteration in time.

We used the finite difference method [22, 23] to solve the transmural pressure equation (2.18) and governing equation (2.39). Specifically, the arterial wall was discretized into segments. Next, a node was placed between two consecutive segments. More than 100 nodes were distributed over the wall length. Then incremental numerical calculations were applied at each node.

For the differential terms included in the transmural pressure equation (2.18) and governing Equations (2.39), forward and central schemes [22] were used to solve time- and space-dependent terms, respectively. Also, the trapezoidal rule [23] was utilized for solving the integral term equation (2.19). Numerical accuracy was first assured by systematically increasing the number of nodes until the results became invariant with further increase. Then the numerical code was validated by comparing the predicted results for steady-state cases with available results from other studies that used different solution methods (Djordjevic and Vukobratovic [13] and Luo and Pedley [10]). Our predictions were qualitatively in agreement with those studies.

The solution procedure is initialized by first defining the shape of the stenotic arterial wall with a small slope in fulfillment of the lubrication theory. Next, with the imposed boundary conditions and numerical schemes in place, the transmural pressure is solved first at all the nodes of the wall. The matrix containing the transmural pressure at each node is inserted into the governing equation (2.39). Then the latter equation is solved for all the nodes at once using the implicit Thomas algorithm [22]. The results obtained provide the new location distribution of nodes at the current time step. Finally, new configuration of the arterial wall is used to update the initial configuration for the next time step. The procedure continues in time until the input number of time steps is achieved or the solution becomes unstable. Stability is controlled by the Courant number 𝐶, defined as 𝑈_max∗Δ𝑡/Δ𝑥<𝐶, where 𝐶=1 and 𝑈_maxis the maximum flow velocity. All the computations are performed using Matlab 6.5 software.

3. Results

In this section, we present the results and discuss the effect of the control parameters, including inlet velocity 𝑈_𝑜, (or equivalently the Reynolds number, Re), inlet transmural pressure 𝑃_𝑡, and arterial wall thickness ℎ. The other parameters are kept constant, including arterial segment length 𝐿, Elastic modulus𝐸, Poisson ratio 𝜗, fluid kinematic viscosity, 𝜈, initial wall deflection relative to nominal radius, 𝛾, and reference pressure 𝑃_𝑜 which is taken as atmospheric. The input data utilized in this study are presented in Table 1, in which 𝛿_𝑜 is the segment radius. These input values have been selected based on previous work and established data on human coronary artery [24–26]. The segment length 𝐿=100 mm, the nominal diameter 𝛿_𝑜=3 mm, and inlet velocity (and correspondingly, the Reynolds number range) are typical in straight sections of adult human coronary arteries [24, 25]. The Young’s modulus (𝐸) and Poison ratio (𝜗) are elastic properties of the wall and are similarly chosen to be typical of human coronary arteries [24–26]. It should however be noted that the elastic properties may vary depending on the extent of disease. The viscosity is a well-established property of blood of normal hematocrit [25]. The transmural pressure 𝑃_𝑡, wall thickness ℎ, and inlet velocity 𝑈_𝑜 are control parameters which do not necessarily have to be constant in a single artery. For example, the thickness ℎ may vary depending on the type of plaque (focal or distributed) and location of plaque and segment of the artery being considered. Inlet pressure depends on the blood pressure which may be patient specific. Inlet velocity and pressure similarly depend on the location of the coronary artery being considered relative to the aorta. The control parameters could therefore be varied in a parametric analysis as done in this paper to investigate the effect of blood pressure, intima thickness, and flow rate. The ramge of values used for the study are typical of human coronary arteries [24–26].

Table 1 

Input data utilized in this study.

In order to compare our results with the steady-state solution of Djordjevic and Vukobratovic [13], we applied the parabolic profile velocity of (2.10) at the tube inlet. Figure 3 shows the time-dependent cross-sectional average velocity 𝑈(𝑡,𝜉) (2.13) calculated at the inlet and at midlength of segment. We applied oscillatory flow velocity to simulate pulsatile blood flow in a coronary artery. It should be noted that the average velocity varies not only with time but also space as fluid flows through the constricted segment of the channel.

Figure 3 

Average velocity profile of pulsatile flow applied at the inlet and midlength of segment.

3.1. Velocity Profile

Figure 4 shows typical predicted profile of magnitude of velocity (inlet velocity 𝑈_𝑜 = 0.2 m/s, 𝑡=0.1 s) in the tube in which the y-axis represents the radial direction from centerline (𝑦=0) to the wall (e.g., 𝑦=1 at inlet), and 𝑥 the longitudinal direction from inlet (𝑥=0) to the outlet (𝑥=6). Minimum velocity is represented with blue in the color code and maximum with red. The result indicates that at every axial location (𝑥) the velocity varies from a minimum at the wall to maximum at the centerline, in the consensus of a developed parabolic profile in a tube. The variation in velocity is a result of fluid viscosity, which allows the velocity relative to the wall to be zero, while reaching a maximum at the axis of the channel. The largest velocity in the tube occurs at midlength between the inlet and the outlet sections where the tube is most constricted (i.e., minimum cross-sectional area), as expected, in order to satisfy mass conservation. These results while not novel, serve to validate the model as it correctly predicts the expected velocity profile in the tube.

Figure 4 

Typical velocity profile at time 𝑡=0.1 shows that the flow is fully developed throughout the channel.

3.2. Effect of Inlet Velocity

Figure 5(a) shows the predicted transmural pressure distribution along the tube for three inlet velocities 𝑈_𝑜=0.1 m/s (thin blue line), 0.2 m/s (purple continuous line) and 0.3 m/s (red dash line). The Reynolds number ranges from 3.2 to 9.6, well within the lubrication theory regime assumed in the analysis. The other parameters such as transmural pressure at the inlet (𝑃_𝑡𝑜 = 180 Pa), and wall thickness ℎ=0.8 mm, are kept constant. The corresponding wall shear stress distribution is presented in Figure 5(b). In general, the transmural pressure decreases (magnitude increases) from the inlet to the trough of the channel and then drops (magnitude increases) in the downstream tube expansion towards the outlet. The slight oscillation in pressure at the end of the segment may be attributed to artifacts of the numerical scheme in satisfying the imposed boundary condition at the tube ends. Figure 5(a) also shows that 𝑃_𝑡 is highest (magnitude is lowest) and is nearly uniform at the highest inlet velocity, 𝑈_𝑜 = 0.3 m/s. The coupling of Reynolds number (inlet velocity) to transmural pressure 𝑝=𝜀^𝑛𝑃(𝜉,𝑦,𝑡) of (2.2) is expressed directly through (2.37) for the parameters 𝑛and 𝜀. Increasing the Re reduces the magnitude of the transmural pressure. Corresponding to the observed intramural pressure trend, Figure 5(b) shows that the wall shear stress is nearly uniform at the highest inlet velocity when the large flow momentum is sufficient to overcome the effect of tube constriction.

(a)

(b)

Figure 5 

Intramural pressure distribution (a), and wall shear stress distribution (b) for inlet velocity 𝐔𝐨 = 0.1 m/s, 0.2 m/s, 0.3 m/s.

Note that all the results presented in Figure 5 (as well as subsequent ones) are characterized by asymmetry between the tube inlet and outlet due to the effect of fluid flow. This observation is consistent with results obtained in previous studies [10, 27].

Figure 6 presents the transient wall deflection corresponding to the pressure and wall shear stress results of Figure 5 above, for inlet velocity 𝑈_𝑜 ranging from 0.1 to 0.3 m/s. In each of the plots, the direction of increasing time is from bottom to top. The results indicate that the amplitude of wall deflection decreases as Re increases. In other words, the wall progressively flattens out as Re increases, in the consensus of the previously observed flattening of the pressure curve (due to significant reduction in magnitude of transmural pressure) and wall shear stress curve at high Re in Figure 5. Figure 6 also shows the wall deflection is the highest towards the half length of the segment. This trend is consistent with the steady-state results of Djordjevic and Vukobratovic [13] and Luo and Pedley [10, 27]. There is significant wall oscillation, especially at the smallest Re. This wall oscillation, when considered in conjuction with the magnitude of the transmural pressure presented in Figure 5(a) suggests possible approach to breakdown in equilibrium.

(a)

(b)

(c)

Figure 6 

Wall shape for inlet velocity 𝑈𝑜 of 0.1 m/s (a), 0.2 m/s (b), 0.3 m/s (c).

3.3. Effect of Inlet Transmural Pressure

Figure 7(a) shows the predicted transmural pressure 𝑃_𝑡 at 𝑡=0.1 s along the segment length for three inlet transmural pressures (𝑃_𝑡𝑜 =180 Pa, 200 Pa, and 220 Pa). The Reynolds number Re=6.4 and wall thickness ℎ=0.8 mm. The corresponding wall shear stress distribution is presented in Figure 7(b). The results indicate that increasing 𝑃_𝑡𝑜 reduces both the magnitude and longitudinal gradient of the transmural pressure along the tube. This finding is consistent with a previous steady state study [27] indicating that high 𝑃_𝑡𝑜 effectively results in a reduction of longitudinal tension in the membrane, or increase in the internal pressure. The net result is a flattening of the membrane as would be described below. The flattened wall in turn leads to a reduction in wall frication, translating to a reduction in wall shear stress as observed in Figure 7(b).

(a)

(b)

Figure 7 

Distribution of pressure (a) and wall shear stress (b), for inlet transmural pressure of 180 Pa, 200 Pa and 220 Pa.

The predicted transient wall deflections corresponding to the results of Figure 7 are presented in Figure 8, with (a) representing 𝑃_𝑡𝑜 = 180 Pa, (b) 200 Pa and (c) 220 Pa. Increasing 𝑃_𝑡𝑜 results in outward movement of membrane (i.e., increased wall deflection). This trend is consistent with the increase in internal pressure resulting from increased 𝑃_𝑡𝑜 and has also been observed in a previous steady-state study [13]. A bulge is predicted at both the upstream and downstream ends of the segment. The former is consistent with that observed in a previous study utilizing higher initial wall constriction [27]. The second downstream bulge, observed here, is not surprising in such a situation of high transmural pressures coupled with small wall constriction. Specifically, the small initial wall deflection implies some degree of symmetry at both ends of the segment.

(a)

(b)

(c)

Figure 8 

Wall shape for inlet transmural pressure, 𝑃𝑡𝑜 = 180 Pa (a), 200 Pa (b), 220 Pa (c).

3.4. Effect of Wall Thickness, ℎ

Figure 9(a) shows the predicted transmural pressure, 𝑃_𝑡, at 0.1 s along the tube for three wall thicknesses ℎ=0.7 mm, 0.8 mm, and 0.9 mm. The inlet parameters are kept constant at Re = 6.4, and 𝑃_𝑡𝑜 = 180 Pa. The corresponding wall shear stress distribution is presented in Figure 9(b). The results show that increase in ℎ reduces the magnitude of the transmural pressure, 𝑃_𝑡. This decrease in 𝑃_𝑡 increases the wall shear stress, progressively destabilizes the membrane, and may ultimately lead to its collapse as described below.

(a)

(b)

Figure 9 

Transmural pressure distribution (a) and wall shear stress distribution (b) for wall thickness varying from 0.7 mm to 0.9 mm.

Figure 10 presents the predicted transient wall deflections as a function of the wall thickness, ℎ. Increase in the thickness from 0.7 mm to 0.8 mm produces a reduction in wall deflection. This trend has also been predicted in a previous steady-state analysis [13]. The present transient analysis shows that further increase in ℎ from 0.8 mm to 0.9 mm dramatically alters the trend. Specifically, although the deflection is still reduced at the trough, the crests now exhibit larger deflections than those for ℎ=0.7 mm and 0.8 mm. The large crest deflections may be attributed to significant increase in the bending force at large ℎ. Note that the wall thickness ℎ appears in the 4th order bending force term of the integrodifferential equation (2.39). This being the dominant term of the equation implies that any slight change in the bending force due to change in ℎ will have significant effect on the wall configuration as observed in Figure 10.

(a)

(b)

(c)

Figure 10 

Predicted wall shape for wall thickness ℎ=0.7 mm (a), 0.8 mm (b), and 0.9 mm (c).

4. Conclusion

A transient analytical model based on the lubrication theory has been developed to study flow-structure interaction in a two-dimensional constricted distensible tube with compliant wall. Previous analytical studies have focused on the steady state. The long-term objective of the present study is to use the transient model developed to provide a relatively simple means of investigating the essential features of stenosed (diseased) vascular remodeling. Such remodeling occurs due to the elasticity of the walls and the pulsatile fluid flow through the vessel. The resulting wall shear stress pattern and wall dynamics have been implicated in the genesis of vascular diseases such as atherosclerosis, as well as plaque rupture that causes to myocardial infarction (heart attack) and may ultimately lead to death. In the present model, the Karman plate theory in the transient mode is used to obtain expressions for the transmural pressure in order to derive a single integrodifferential equation for the wall behavior. The equation was then solved by the finite difference method. Parametric studies were conducted on the effect of three parameters: Reynolds number, Re (or correspondingly the inlet velocity, 𝑈_𝑜), inlet transmural pressure, 𝑃_𝑡𝑜, and wall thickness, ℎ.

The primary findings of the analytical study for the specific conditions considered may be summarized thus. (i)Increase in Reynolds number while keeping 𝑃_𝑡𝑜 and ℎ constant results in a reduction in the magnitude of the transmural pressure. This reduction in 𝑃_𝑡 in turn results in a reduction in both the wall deflection and wall shear stress. The wall is most vulnerable to oscillations at low Re when the other system parameters (𝑃_𝑡𝑜,ℎ) are fixed. The wall deflections are also generally larger upstream of the trough (minimum diameter at midlength) than downstream. (ii)Increase in the inlet transmural pressure 𝑃_𝑡𝑜 while keeping ℎ and Re constant causes a reduction in longitudinal tension and a flattening of the wall. This trend results in a reduction in the wall friction and wall shear stress. (iii)Increase in wall thickness ℎ while keeping 𝑃_𝑡𝑜 and Re constant reduces the wall deflections for ℎ below a threshold (here ℎ=0.8 mm). Increasing ℎ beyond this threshold leads to large deflection especially at the ends and possible breakdown of the wall. This collapse of the membrane is due to large bending force that develops at the ends at large wall thickness.

The results at specific time duration seem to agree generally with previous steady-state results [10, 12, 13]. They are also generally in the consensus of medical studies on cardiovascular diseases. For example, the low shear stress resulting from increase in the inlet transmural pressure 𝑃_𝑡𝑜 implies that the latter may indeed promote atherogenesis. The latter has been found to be focal in coronary arteries and located in regions of low wall shear stress. These are also regions of platelet accumulation. The observed deflection is larger upstream than downstream of the trough. This trend is consistent with observations that plaques typically rupture at upstream shoulders, that is, upstream of the region of minimum vascular diameter. Increase in wall thickness may be directly related to intimal thickening due to vascular inflammation. The present results indicate that such thickening may lead to possible breakdown or rupture of the wall and plaque especially near the ends. Plaque rupture causes thrombosis and may lead to heart attack and death.

Our analysis has demonstrated that with proper assumptions the transient model can be used to study the unsteady flow-structure interaction phenomena in a constricted compliant channel. The results obtained from the investigation are consistent with previous analytical studies [10, 12, 13]. The study is currently limited to system behavior in the relatively low Reynolds number range 3.2 to 9.6 in order to simulate flow in the coronary arteries that are most vulnerable to vascular plaque growth and rupture. Subsequent efforts will be extended to higher Reynolds numbers in order to investigate other interesting phenomena such as self-excited oscillations for direct comparison with available experimental data that were conveniently conducted at high Re. In addition, it is desirable to improve the model by extending to larger initial membrane curvature.

Next, it would be interesting to perform stability analysis of the improved model to investigate membrane self-excited oscillation [10, 27, 28]. The investigation will help to optimize the model by determining the role each controlling parameter plays in order to focus on the most important one in analyzing the fluid-structure phenomena. Ultimately, the model would be extended to investigate the fully coupled three-dimensional model. The investigation of the fully three-dimensional geometry [10, 27] in order to investigate the subtleties compared to the 2D model considered in the present study. Such subtleties include snap-through buckling of the tube wall, flow division into multiple lobes which remains open during the buckling, and buckling with higher circumferential wave numbers for increased upstream pressure.