Abstract
The accumulation of microscale materials at solid–fluid interfaces in biological channels is often the initial stage of certain growth processes, which are present in some forms of atherosclerosis. The objective of this work is to develop a relatively simple model for such accumulation, which researchers can use to qualitatively guide their analyses. Specifically, the approach is to construct rate equations for the accumulation at the solid–fluid interface as a function of the intensity of the shear stress. The accumulation of material subsequently reduces the crosssectional area of the channel until the fluidinduced shear stress at the solid–fluid interface reaches a critical value, which terminates the accumulation rate. Characteristics of the model are explored analytically and numerically.
1. Introduction
The primary objective of this work is to model a ‘generic’ solid–fluid interface accumulation in biological channels, which contains flowing fluids with suspensions, and to develop a model that is relatively easy to evaulate. In some forms of growth, the first stage is attributed to accumulation of microscale suspensions in the fluid which adhere to the flow boundaries (walls). It is this stage, accumulation, that is the focus of this work. For example, one important application that motivates the present analysis is plaque growth because of high lowdensity lipoproteins (LDL) content. The phenomena of plaque buildup are thought to be due to a relatively high concentration of microscale suspensions (LDL particles) in blood. Plaques with high risk of rupture are termed vulnerable. Atherosclerotic plaque formation involves: (I) adhesion of monocytes (essentially larger suspensions) to the endothelial surface, which is controlled by the adhesion molecules stimulated by the excess LDL, as well as the oxygen content and the intensity of the blood flow, (II) penetration of the monocytes into the intima and subsequent tissue inflammation and (III) rupture of the plaque, accompanied by some level of thrombus formation and possible subsequent occlusive thrombosis. For surveys of plaquerelated work, see the earlier studies [1–11]. The mechanisms involved in the early part of plaque formation (stage (I)) have not been extensively studied, although some qualitative studies have been carried out (e.g. [12,13]). Another related problem that motivates the upcoming analysis is calcification in the aortic valve. The deposition of calcium and fatty deposits on the valve, in relation to the shear stress on the leaflet wall, and how this leads to different tendencies of growth on the sides of the leaflet (facing the aorta or left ventricles) is an open question of increasing interest (e.g. [14–17]).
The approach in the present work is to develop rate equations for the accumulation of suspended microscale material at the solid–fluid interface. In the model, the intensity of the shear stress at the interface dictates the rate of accumulation. As a channel cross section narrows, the flow rate becomes more intense, thus naturally limiting the accumulation. Both analytical and numerical approaches are undertaken to determine the model's characteristics.
2. The basic model
Consider an idealized channel with a circular cross section of (initial) area with an initial solid–fluid radius R_{0} (figure 1). The objective is to describe the mechanism by which A_{0} changes, via a reduction in R_{0}, owing to wall accumulation caused by microscale deposits building up onto the channel's interior surface. The interior channel (solid–fluid interface) radius, denoted R, changes over time (R(t = 0) = R_{0}). As a simplification, we assume that, at a given location along the channel, the accumulation is radially symmetric, and that at a given longitudinal (z) location, at the wall, r = R, and that we have a velocity profile given by a classical channelflow of the form 2.1where v_{max} is the centreline velocity. For fully developed laminar flow, q = 2, while for increasing q one characterizes, phenomenologically, progressively turbulent flow (q ≥ 2). The shear stress is given by 2.2The (reaction) stress at the wall (at r = R) is . We will further assume that the overall flow rate is assumed constant 2.3
One can show that 2.4
The stress at the wall becomes 2.5
We have the following observations:

— increasing µ, Q_{0} or q increases the stress at the wall (τ_{w}),

— increasing q leads to an increasingly more blunted flow profile and

— decreasing R increases the stress at the wall (τ_{w}).
Remark.
In the next few sections, we will assume that the flow profile exponent (q) is independent of velocity. However, later in the presentation, we will relax this restriction and correlate q to the centreline Reynolds’ number (Re).
3. Rate of interface accumulation
It is assumed that the tendency for material to adhere to the wall is controlled by the intensity of the shear stress near the wall. Essentially, higher shear stresses reduce the likelihood of material adhering to the wall, whereas lower shear stresses increase the tendency of material to adhere to the wall. Accordingly, we consider the accumulation rate to be proportional to the nondimensionalized fractional difference between the shear stress at the wall (at a fixed zlocation) and the critical ‘detachment’ stress (τ* > 0). Specifically, for the reduction of the solid–fluid radius 3.1where R_{0} – R(t) is the reduction in radius (figure 1) and η is a rate constant representing the accumulation per unit time. As the velocity increases (increasing the shear stress), the microscale material in the fluid is less likely to adhere. We have the following observations:

— increasing τ* increases the rate of adhesion,

— the unilateral limiter τ_{w} ≥ τ* in equation (3.1) shuts off the accumulation at the solid–fluid interface,

— the rate of accumulation decreases with increasing µ, Q_{0} and q and

— the rate of accumulation decreases with decreasing R.
Remark 3.1.
Assuming that τ_{w} < τ*, the accumulation equation has the following form: 3.2where and a_{2} = η. If we linearize ℱ(R) about R_{0} = R(t = 0), and we obtain 3.3where the superscript L indicates a linearized value and and . Using standard techniques (superposing homogeneous and particular solutions), we obtain 3.4Clearly, the rate of decay (of R) is dictated by η, and that decay is expected to be exponential. Beyond that, the linearized solution provides little insight into the character of general model, which is treated numerically next.
Remark 3.2.
The steadystate value of the solid–fluid interface, denoted R^{ss}, can be determined by setting dR/dt = 0, leaving 3.5which illustrates

— the steadystate value R^{ss} is independent of η,

— increasing the detachment stress threshold (τ*) leads to more accumulation (reduction of R^{ss}), with an inverse cubic dependency and

— increasing μ, Q_{0} and q leads to less accumulation (larger R^{ss}) with a cubic dependency.
Remark 3.3.
One can compute the amount of change of R^{ss} owing to changes in the parameters from the sensitivities to the parameters in the system and Furthermore, the ratios of the sensitivities provide information on the relative sensitivity of one parameter to another, for example and
This allows one to determine the amount of change of one parameter, say δQ_{0} that would be needed to make a comparable change δR^{ss} due to say δμ by computing 3.6The procedure is virtually identical for the other parameters in the system.
4. Direct timetransient numerical simulation
The general form of the equation 4.1is solved using an explicit forward Euler time integration of the form 4.2The simulation time was set to 10 years.^{1} An extremely small (relative to the total simulation time) timestep size of 5 × 10^{–4} ϕ, where ϕ = 3600 × 24 × 365 is the number of seconds in a year, was used. Further reductions of the timestep size produced no notable changes in the results, thus the solutions generated can be considered to have negligible numerical error. The following physical parameters were used:

— R(t = 0) = R_{0} = 0.01 m,

— m^{3} s^{−1}, where v_{m}(t = 0) = 0.1 m s^{−1} in the mean velocity,

— μ = 0.003 Pa s^{−1},

— τ* = 1 kPa,

— ρ = 10^{3} kg m^{−3} and

—
The results are shown in figure 2. As indicated in figure 2, eventually, the accumulation slows, then terminates, once the channel narrows sufficiently to raise the fluidinduced shear stress to exceed the threshold value of τ*. For the set of parameters chosen, the initial centreline velocity is v_{max} = 0.2 m s^{−1} and ramps up to v_{max} = 82.207 m s^{−1} owing to the reduction of the cross section and the constant volumetric flow rate Q = Q_{0} (figure 3). Realistically, as the crosssectional narrows, the flow becomes relatively more turbulent, i.e. the Reynolds’ number increases. The exponent q must then also change from a quadratic profile a more blunted profile. To incorporate this effect, we introduce a dependency of q on the centreline Reynolds number shortly. This will have the effect of further accumulation rate limitation. Clearly, for increasing η, the rate of accumulation increases, and the time to steadystate decreases. While figure 2 illustrates the basic behaviour with respect to variations in the accumulation rate parameter, there are clearly other possible parameters to vary, and we refer the reader to §3 pertaining to a sensitivity analysis.
Remark.
We note that one can write the time needed to reach a given value of R(t*) = R* in terms of the following integral: 4.3which can be numerically integrated, however, this timescale can also be inferred easily from figure 2.
5. Model extensions: rate limiting effects—blunted velocity profile evolution
Clearly as the Reynolds number increases, the velocity profile will change from a quadratic (q = 2) to a more blunted profile (), which represents, phenomenologically, turbulent (inertiadominated) behaviour (figure 4).
5.1. Incorporating profile changes
The effect of a changing profile is described by representing q by a linear function of the centreline Reynolds’ number (Re) 5.1where and c_{1} and c_{2} are constants. Models of this type, linking the the profile exponent (q) to the centreline Reynolds’ number (Re), are quite wellestablished (e.g. [18]). Usually, and c_{2} ≈ 2, and in the limit we have, for c_{1} = 0 and c_{2} = 2, laminar flow (q = 2). For the general case, combining equation (2.4) with equation (5.1) and the definition of the centreline Reynolds’ number, we obtain a quadratic relationship for q, 5.2where . This quadratic relationship can be solved in closed form for q to yield^{2} 5.3Clearly, the larger root is the physically correct choice. Q(Re) is clearly a function of R^{–1} and decreasing R increases q.
5.2. Numerical examples
The general form of the equation 5.4is solved using a forward Euler time integration, . The same parameters as for the q ≠ q(Re) case were used, with the two additional parameters for the change of flow profile with Reynolds' number, c_{1} = 10^{–2} and c_{2} = 2. Figures 5 and 6 illustrate the rate limiting effect of the blunting of the flow profile (increasingly more turbulent flow), relative to the case where q is not a function of the flow velocity (figures 2 and 3). For the set of parameters chosen, the initial centreline velocity is v_{max} = 0.2 m s^{−1} and ramps up to v_{max} = 9.716 m s^{−1}, which is much more realistic than the constant q model which predicts v_{max} = 82.207 m s^{−1}. Numerically, the solution of this more relatively complex model q = q(Re) requires virtually no additional effort relative to q ≠ q(Re). One could attempt to analytically extract the steadystate value of R by setting and combining it with equation (5.3), leaving 5.5
However, this leads to a highly nonlinear equation, with multiple roots, that would have to be solved numerically. Clearly, it is easier to simply track the steadystate value of R from figure 5.
6. Conclusion
As indicated at the beginning of this paper, the primary objective of this communication was to investigate solid–fluid interface accumulation on the walls of biological channels and to develop a model that researchers in the field can easily implement and use to guide experiments. Semianalytical relationships were also developed which characterize accumulation and flowinduced rate limitations, owing to the more intense flows for smaller (occluded) cross sections. This is useful because of the longterm character of experiments involved with tracking accumulation. Future work should concentrate on the development of a subsequent growth model, involving both adhesion as well as constitutive growth (e.g. [19–23]). It is also important to explore how growth depends on temporal/spatial gradients of wall shear stress. These processes probably involve strongly coupled diffusive, chemical effects and thermal effects, and we refer the reader to Markenscoff [24–26] for in depth mathematical analysis of such coupled systems. There are conflicting reports on whether high or low shear stress correlates with growth. Certainly, there may be scenarios where growth is controlled by low shear stress in one phase and high shear stress in another. Clearly, the mechanism for biological attachment is quite involved, and the simple notion used in this paper of a critical attachment/detachment stress threshold may be inadequate. Therefore, models building upon results such as those found in Hermanowicz [27–29], Sawyer & Hermanowicz [30] and Yoon & Mofrad [31] may prove quite useful in this regard.
The presented analysis and model can provide a useful guide to designing and interpreting experiments, which can take years. However, while the model can provide qualitative a priori information for further computationally intensive largescale simulations, extensions are invariably going to require complex spatial discretization of the system under analysis and could also entail resolving particle–fluid interaction. The number of research areas involving particles in a fluid undergoing various coupled processes is immense, and it would be futile to attempt to catalogue all of the various applications. However, a common characteristic of such systems is that the various physical fields (thermal, mechanical, chemical, etc.) are strongly coupled, with particles that tend to agglomerate (cluster). In Zohdi [32], a flexible and robust solution strategy was developed to resolve coupled systems comprising large groups of flowing particles embedded within a fluid, based on agglomeration models found in Zohdi [33]. In that analysis, particles were surrounded by a continuous interstitial fluid which is assumed to obey the compressible Navier–Stokes equations. Thermal effects were also considered. Such particle/fluid systems are strongly coupled because of the mechanical forces induced by the fluid onto the particles and viceversa. Because the coupling of the various particle and fluid fields can dramatically change over the course of a flow process, a primary focus of that work was the development of a recursive ‘staggering’ solution scheme, whereby the timesteps were adaptively adjusted to control the error associated with the incomplete resolution of the coupled interaction between the various solid particulate and continuum fluid fields. The approach is straightforward and can be easily incorporated within any standard computational fluid mechanics code based on finite difference, finiteelement, finite volume or discreteelement discretization, for example, those developed in earlier studies [34–39].
Footnotes
 Received October 8, 2013.
 Accepted November 7, 2013.
 © 2013 The Author(s) Published by the Royal Society. All rights reserved.