The effect of cartilaginous rings on particle deposition by convection and Brownian diffusion

The deposition of spherical nanoparticles by convection and Brownian diffusion in a pipe with a cartilaginous ring structure is studied. Analytical results for a fully developed flow are found for small amplitude rings using the interactive boundary layer theory. It is found that the local deposition rate is at maximum at a position approximately one twelfth of the spacing between the rings before the minimum cross section of the tube. For larger ring amplitudes the problem is solved numerically and separation then takes place in the depressions between the rings, and maximum deposition is found at the point of reattachment of the flow approximately at the same point as in the analytical theory. Cumulative deposition results are also provided with larger deposition rates with the inclusion of the cartilaginous rings. Deposition results for a developing flow are also provided. For the same volume flux as for fully developed flow the deposition is about 25% larger. In general conclusions about the position of maximum deposition rate from the analytic theory of fully developed flow also applies qualitatively to the case of developing flow.


INTRODUCTION
Usage of nanoparticles in material design enables the development of new superior products as well as improvement of existing ones, but these tiny particles are potentially toxic and pose substantial health risks [1].
There is an increasing literature that discusses the ability of nanoparticles to cause adverse effects in cellular functions [2].Depending on the field of application, nanoparticles may enter the body via inhalation, dermal, oral or injection exposure routes.Carbon nanotubes, for instance, have a diameter of approximately 10 nm with a very much larger length of 1-10 µm, and if inhaled then, they may probably quite easily move through the airways down to the alveoli.Once in the alveoli, they can get stuck and initiate different diseases, in a similar manner as asbestos fibres.Thus, knowledge of transport and deposition properties of aerosol particles in lung flows is essential.Such knowledge is also useful in the optimization of drug delivery with pharmaceutical aerosols.
Of special interest is the effect of cartilaginous rings, located in the upper generations of the human airways, on nanoparticle deposition caused by convection and diffusion.These rings give the airway wall an uneven surface, with ring forms protruding into the airway lumen, and with depressions between the rings.These regular and symmetric structures are expected to influence turbulent as well as laminar airflow directly passing them, which could potentially affect deposition in this region.Measurements of a physical model with cartilaginous rings have been done by Zhang and Finlay [3] using particles of size about 5 µm.They find an increase in deposition of about 20% due to cartilaginous rings.
In the present study, we examine the diffusion of smaller spherical particles and their deposition on the walls of an airway with cartilaginous rings, while leaving the corresponding problem for fibres for future studies.The results provided by the present analysis for spherical particles can then be considered as an upper limit for the deposition of fibres.
One method to determine particle deposition in the diffusion limit is to solve the Langevin equation with a stochastic force on the particle in a given flow field cal-culated from the Navier-Stokes equations.To find deposition statistics, however, a large number of particles need to be launched [4].An alternative approach is to consider the equation describing the evolution of the probability density, the Fokker-Planck equation, which has the advantage of directly providing the statistical averages [5].The disadvantage of this method is that for large Peclet-numbers (Pe), implying small diffusion, very thin boundary layers develop, and these are difficult to resolve numerically.The Fokker-Planck equation has the same form as the convective-diffusion equation, describing the evolution of the concentration of particles.In the present paper, the particle deposition is derived by combining the Navier-Stokes equations and the convective-diffusion equation.This methodology has previously been applied by Ingham [6][7][8] and Martonen et al. [9] for deposition in smooth-walled tubes.
A lot of research has been done using a similar set of equations describing in heat transfer and mass transfer with wavy walls that often occurs in engineering applications.In these applications however the Pe-number is often much smaller, of the same order as the Re-number.
Here we consider Pe much larger than the Re-number.We also believe that these general results are not so well known in connection to the biological application considered in the present paper.
The paper is organized as follows.In Section 2, an approximate analytic solution of the Navier-Stokes equation is presented together with solutions of the convective-diffusion equation.To solve the problem for the flow with a weakly perturbed sinusoidal boundary, a linearized interactive boundary layer theory is employed [10].To start with, and as a zero-order solution, a fully developed pipe flow is chosen.The analytic solution has the advantages that it provides a better understanding of the flow and that it introduces the relevant scaling of the problem.Furthermore, it can be used to verify the corresponding numerical results, to be presented in Subsection 3.1 Due to the previously mentioned numerical difficulties with very large Pe, we only consider light breathing conditions, which implies large but modest Pe.In Subsection 3.2, a developing flow is considered.For this case, it turns out to be handy to take a numerical approach, since an analytical treatment is difficult, even without rings.To validate the numerical work, the results are finally compared with approximate treatments of entry flow in smooth-walled tubes presented by Ingham [7][8] and Martonen et al. [9].Here two rather different configurations of cartilaginous rings are considered.
In a real respiratory system neither fully developed nor developing axisymmetric tube occurs.It is rather a combination of the two since after an airway bifurcation approximately half the flow attaches a new wall creating a new developing thin velocity boundary layer.However we believe that the two extreme cases considered in this paper gives a first estimate of the number of deposited particles in the real case.

The Interactive Boundary Layer Theory
The axisymmetric geometry in Figure 1 is considered, for which the flow may be described by the dimensionless Navier-Stokes equations according to Here u is the stream-wise velocity and v the radial velocity, and p is the pressure.In addition, the convectivediffusion equation is introduced in the form in order to study the diffusion of particles embedded in the fluid.In (2), c can be considered as concentration of particles, while in a Fokker-Planck context as a probability density function.
The Reynolds number Re and the Peclet number Pe are now defined as where 0 U is the maximum velocity for the fully devel- oped pipe flow, U is the mean velocity at the inlet for the developing flow,is the viscosity, D is the diffusion coefficient, and a is the inlet pipe radius., where λ is the mean free-path and L is a characteristic length scale of the problem.For small Kn the Navier-Stokes equations can be derived from an expansion in small Kn-numbers of the more fundamental gas kinetic theory (Chapman and Cowling, [11]).For the flow in a pipe the radius can be considered as a characteristic length and then the Knudsen number is very small in the present application.For the flow past a spherical particle with radius of the order of nano-meters, the Knudsen number is however not small and the Navier-Stokes equations are not valid.In a layer close to the particle, of the order of the mean-free path of the air-molecules, the more fundamental gas kinetic theory must be applied.The role of the Navier-Stokes equations is then to provide with "outer" local values of the velocity field for the free-molecular flow around the particle.The drag force and the drag coefficient of the spherical particle are only known theoretically for small Kn numbers (continuum limit) and large Kn numbers (freemolecule limit).For the continuum limit the drag coefficient is where f d is the particle diameter.For the intermediate range of Kn numbers an empirical correction, known as the Cunningham factor is introduced as 1 (2.34 1.05exp( 0.39 )) The expression for the drag coefficient is then which provides a value for the drag coefficient for all Knudsen numbers, Kn / f d   .For Brownian diffusion the diffusion coefficient D is then related to the drag as where is Boltzmanns constant and T is the temperature.The wavy boundary is described in terms of the dimensionless amplitude , defined as * a    .
In order to analyze the flow problem analytically, an interactive boundary layer theory is applied to the fluid flow problem.Following the presentation by Sobey [10], the solution of the Navier-Stokes equation is divided into two parts, the outer part describing the core flow and the inner one describing the flow in the neighbourhood of the boundary.The present formulation of the theory is applicable when the typical length scale of the wall perturbation is of the same order as the pipe radius.In this so called "fine indentation limit" a balancing of the order of magnitude of terms gives that the wave amplitude scales with the Reynolds number as The core flow is then described by the small  and large re-number expansions in which the leading order terms correspond to unperturbed pipe flow, and the second terms are corrections due to the finite amplitude of the wavy boundary.In the correction terms note the relationship between the amplitude and the Reynolds number given in (8).For a description of the boundary layer, a coordinate R is introduced in the following manner: within the present formulation, the velocities in the boundary layer are scaled as .
The boundary layer equation for the axial velocity component can be defined as where the pressure gradient is only a function of the stream-wise coordinate and the pressure scales as It is convenient to rewrite (12) in terms of the streamfunction  of the boundary layer as Another wall coordinate is introduced as is the shape of the wall and is a parameter controlling the amplitude of the wavy wall.The boundary layer (12) then becomes where the asymptotic matching condition with the outer core solution is given by and the boundary condition at the wall To proceed analytically, the boundary layer ( 17) is solved for the case of a small perturbation of the cylindrical wall so that << 1.
This means that solutions corresponding to flow separation cannot be strictly described, however increasing from small values give a first indication of where the flow may separate.For larger wall perturbations   (17) has to be solved numerically together with boundary conditions and matching with the core solution and the upstream solution.For small µ we assume that the solution of the boundary layer equation can be written as 2 ( , ).
The matching with the outer core solution is then expressed as

( ) . f x asY Y
After linearization the boundary layer, (17) becomes which can be solved by introducing the Fourier transform ˆexp( ) ikx dx A solution for the first derivative is of the form where ˆ( ) W k is a general function of k and Ai is the Airy-function.Matching with the core solution requires that ˆ2 ( ) which determines the function ˆ( ) W k and we find The boundary layer velocity close to the boundary, therefore, is given by the Taylor expansion After inverting the Fourier transform, the velocity becomes For a specific form of the boundary f(x), the function g(x) in (28) can be calculated.An interpretation of this term is that it corresponds to the wall shear stress introduced by the wavy boundary.The minimum of this function is therefore given to the position of minimum wall shear stress.Thus, in a situation in which µ increases from zero to small but finite values, this minimum gives an indication of where the flow first separates.The corresponding wall normal velocity to (28) is The expanded results (28-29) are of importance for the large Pe solution of the convective-diffusion equation since, as will be shown below, only the velocity components close to the boundary are required.

The Convective-Diffusion Equation
Next, consider the convective-diffusion equation Since we consider the limit of large Pe, the lowest order approximation to (30) is simply This suggests that the concentration is preserved along streamlines.For absorption of particles reaching the boundary, however, the required boundary condition is 0 c  , which cannot be fulfilled by solutions of (31), therefore, a boundary layer analysis of (30) is required.Using the boundary variables introduced in the previous section, the diffusion equation can be written in the form 773 where the second derivative of c with respect to x is neglected as is a common praxis in boundary layer analysis.
For the application in focus, is usually very small, so we consider this to be the case for the rest of this paper.Introducing a boundary layer variable By Taylor-expanding the arguments of the velocity components and balancing the order of magnitude of the terms it is found that 1/3 implying that only the approximate form of the velocities in (28-29) are required.
A similarity solution exists with a similarity variable of the type and with a corresponding differen- Choosing the quantity inside the bracket to be 3, the equation for the concentration boundary layer thickness h(x) is obtained as A solution corresponding to zero concentration boundary layer thickness ( ) h x at entry x = 0 is then given by For µ = 0 we recover the result for convective-diffusion in ordinary pipes.The equation for the concentration, (34), then becomes The solution can be written in terms of the incomplete gamma function with solution fulfilling the boundary condition Of special interest is the particle deposition due to diffusion, which corresponds to the normal component of the particle diffusion flux at the wall.In non-dimensional form this is Neglecting terms of order 2 ( ) O  we find the following expression for the deposition: As an example of a possible cartilaginous wall-shape, we consider the periodic function Here, the parameter k 0 is measure of the spacing between the rings.A convenient lowest order approximation for the analytical analysis is the first two Fourier components of (41), i.e., inspection of bronchoscopic videos shows that for the application of cartilaginous rings, these simple wall shapes are quite reasonable, at least in some parts of the tracheobronchial tree and the shape given by ( 41) is similar to that studied by Martonen et al. [12].The differences between the shape (42) and the shape considered by Martonen will be discussed in Subsection 3.2.For x > 0 the inversion of the Fourier transforms, needed for the calculation of the velocity components, can be found by deforming the complex k-integration contour into the complex plane around a branch-cut chosen along the positive imaginary k-axis.For x < 0 the integration contour can be closed in the lower half complex k-plane yielding that it is zero.This means that there is no perturbation of the boundary layer flow upstream of x = 0.For x > 0 the following result is obtained after inversion of the Fourier integrals where g(x) is defined by is a Lommel function defined in Abramowitz and Stegun [12].For large x the first term in g(x) dominates and is represented by a simple sinusoidal form.
In Figure 2 the function g(x) is plotted and compared with the shape of the boundary f(x).It is noted that both functions have the same period with g(x) ahead of f(x) with a phase that for large x becomes .Since g(x) is the correction of the wall shear stress due to the wavy wall, this result can be used to pinpoint the position where separation first occurs.This position is where g(x) has a minimum, i.e., before the minimum of f(x).
The behaviour of g(x) is also essential for the particle deposition rate since, as will be shown below, g(x) is at its maximum where the concentration boundary layer thickness is at its minimum.Correspondingly, the concentration boundary layer thickness is at its maximum at the point of first separation.From (36) the general expression for the concentration boundary layer thickness can be derived as where it can be noted that the expansion ( 46) is invalid for very small x.The local deposition of particles can be obtained from the normal diffusive flux at the wall, given by 0 Now, by inserting (45) into (47) the local flux is given by or with numerical values Please observe that the scaling with Pe is the same as for ordinary pipe flow ( Since the largest deposition occurs for the minimum thickness of h(x), we note from (44-45) that the largest deposition is localized to the region just before the maximum of f(x), which for the wall geometry defined in (42) is in fact approximately ahead of the maximum of f(x).It follows that the minimum deposition occurs at the point of first separation.
Using expression (44), an analytic expression for the integral in (45) can be obtained as Combining this with (46) the local deposition rate given by (48) can be calculated.In Figure 3, the local deposition rate is plotted as a function of x for a finite and a zero value of .From this figure it is clearly seen that the maximum local deposition is slightly ahead of the minimum radius of the pipe.
It is also of interest to calculate the cumulative depo- sition after a length X by performing the surface integral over the azimuth angle and the x-coordinate.This can be done analytically, but the expression is rather complicated and is omitted for brevity.The cumulative deposition after length X is instead written as To get an expression for the probability of a particle to adhere to the boundary after a length X, (48) is divided by the inlet flux, which in dimensionless quantities is equal to 2 so For cylinder tube flow  a rather neat expression This result is valid for a cylindrical smooth-walled tube, and is in agreement with the result given by Ingham [6] for a fully developed parabolic velocity profile.Analysis of ( 49) and (50) shows that the total deposition for finite  is somewhat greater than for the case of an ordinary pipe, and effects from the waviness of the wall are clearly revealed, as exemplified in Figure 4 where the cumulative deposition given by ( 50) is plotted for different amplitudes of the wall perturbation.The values of the dimensionless numbers are Re = 450 and Pe = 100000, corresponding approximately to light breathing conditions in generation 4 [12] with a particle diameter of 10 nm and pipe radius 2.25 mm.Results are shown for an ordinary smooth pipe, and for ring configurations with small ring-amplitudes  and  where we note the relation with  as Re  .The spacing of the rings is chosen such that there are five rings within a length of 6 radii, corresponding to the wavenumber 0 5 k  .

Comparison of Theory with Numerical Solution for Parabolic Flow
The approximate theoretical analysis in Section 2 can only provide results for small ring-amplitudes without flow separation.For larger amplitudes it is necessary to solve the equations numerically.Hence, a numerical solution of (1-2) is considered.The commercial software Comsol Multiphysics 3.5 is used to serve this purpose, since this code enables the Navier-Stokes and convective-diffusion equations to be solved simultaneously.In Comsol, the finite element method is used, and an adaptive mesh-refinement method is applied.For the implementation of the geometry and for post-processing, Comsol script or Matlab are used, both of which can be conveniently applied together with Comsol Multiphysics.To start with, the equations for the wavy wall are solved (42) and then the theory presented in the previous section is used to validate the numerical results.First, the theoretical (51) and numerical results for smooth-walled pipe flow are compared.The agreement is very good for Pe that equals to 100000 and with a Re of 450, despite the boundary layer approximation used in the theoretical analysis, see Figure 5. Next, the rings with a small amplitude of 0.025diameters is studied.
The agreement between the local deposition rate as calculated from the numerical method and the analytic method presented in the former section is rather good, see Figure 6.The agreement between the methods for the corresponding cumulative deposition is also relatively good; where the deviation is most notable for small and large X, see Figure 7.It is therefore concluded that the numerical results are reliable and we proceed with considering larger values of the ring-amplitudes.
Consider a ring-amplitude of .1,mimicking the size of amplitudes used in other studies of effects of the cartilaginous rings, e.g.(Martonen et al. 1994) [12].There is a qualitative difference in the flow pattern between this case and that for small .For he flow separates in the region between the rings, see the lower plot in Figure 8, while for , there is no separation, see upper plot in Figure 8.
A rough criterion for separation can be obtained from the theory developed in Subsection 2.1.From the expression for the velocity profile at the wall g(x) we find for large x that the wall shear stress becomes For the present case this gives the critical amplitude of about 0.06 which is quite close to the estimate from the numerical analysis.
In Figure 8 the red area corresponds to high concentration c while the blue area corresponds to regions of lower concentration.Remarkably, but in accordance with the theory, the deposition is not larger in the separated regions, instead the maximum deposition occurs at the point of reattachment, i.e., just before minimum contraction.This is also in agreement with the theory developed in Subsection 2.2 where the largest deposition is obtained at the maxima of the function g(x) and the lowest deposition at the point of minimum g(x), the latter corresponding to the position of first separation.In Figure 9 the cumulative deposition for amplitudes .1 and are shown.It is noted that the largest increase in deposition compared to a smooth pipe is in the region before the first ring, where the deposition for the case  is about 40% higher than for the smooth-walled tube.

Analytic and Numerical Results for Developing Flow
Since after an airway bifurcation approximately half the flow in each branch will have a developing character, we here consider a developing flow in an axisymmetric tube with cartilaginous rings.For this case an analytic approach is nontrivial, even without rings; hence, the case with cartilaginous rings is only treated numerically.Neglecting the rings, different approximate results have been presented in the literature.Ingham [7,8] applies an integral approach for the entry flow and solves the convective-diffusion equation numerically.An approximate expression for the deposition, similar to (51), is provided in the form  which is valid for small X/R < 0.02.A similar approximate analytical result using a more heuristic method has been given by Martonen et al. [9] as We note the difference in the scaling for X/Re in (55) and (56).It is therefore of interest to consider this problem numerically where the applied conditions at the inlet correspond to uniform flow and the volume flux is the same as for the parabolic flow case with Re = 450 and Pe = 100000.The numerical result obtained fits well in between the results of Ingham and Martonen as shown in Figure 10, where the solid red line corresponds to the numerical solution and the dash-dotted curve corresponds to the result (55) and the dashed line corresponds to the result (56).For comparison the analytical results for the fully developed case is also shown as dotted.It is seen that the deposition is about 25% larger for the developing flow as compared to the parabolic flow.
Since there is a disagreement in the power between (55) and (56), the power for the numerical solution is calculated to 0.5437 for large X, which is more close to the value 5 / 9 0.5556  provided by Ingham [8] than the value 1/2 given by Martonen et al. [9].On the other hand, the results by Martonen et al. [9], obtained from (56), have a good conformity near the entrance of the tube, which is as expected since the flow is of here of Blasius type, which has the power 1/2 (Schlichting et al. [13]).
Next we consider the total deposition for developing flow for rings with amplitudes z = 0.1 and .See Figure 11.The results show that the largest increase in deposition occurs before the first cartilage ring with an increase of the order of 25% for  and 35% for when compared to the case of a smooth pipe.The total deposition is as for the fully developed case also increased considerably by increasing the amplitude from to  Considering other shapes of the cartilaginous rings is also of interest.The shape studied by Martonen et al. [12] is presented in Figure 12.To compare this case with the case studied above we consider the same number of rings.This means that the region between the rings is larger than for the case of the sinusoidal shape discussed above.Since these regions have lower deposition rates than in the region of smaller cross section the overall deposition is smaller than for the corresponding sinusoidal shapes see Figure 13.

CONCLUSIONS AND DISCUSSION
The deposition of nano-sized spherical particles on thewalls of a cylindrical tube with periodically spaced cartilaginous rings has been explored.For fully developed flow and for small ring amplitudes analytic formulas are derived.It is found that the rings increase the deposition rate and it is shown analytically that the largest deposition occurs at a position approximately /12 before the minimum cross section where  is the spacing between the rings.These results are validated with numerical computations.For larger amplitude of the rings  > 0.06 corresponding to light breathing conditions, the flow separates and the largest local deposition is at the point of reattachment.Although it has a strong influence on the variation in deposition along the tube the total deposition is marginally effected by the amplitude of the rings for small amplitudes  < 0.10.For slightly larger values of the amplitude  = 0.05, the total deposition is however increased considerably.
For the case of a developing smooth-wall flow the deposition rates are 25% larger than for fully developed flow while the additional deposition effect caused by the cartilaginous rings is of the same order as for the fully-developed flow.The rate dependency can be of interest for flow in a tube downstream a bifurcation where approximately half of the flow attaches to a new wall.A comparison with the ring shapes considered by Martonen et al. [12] also illustrates another effect, the effect of the ratio of ring radius to the length between the rings.Considering the same number of rings as the sinusoidal shapes, the distance between the rings is larger and therefore from theory we know that these regions of separated flow have lower deposition rates than in the regions close to the rings, and therefore the overall deposition is smaller.
Future work involves other breathing conditions especially heavy breathing conditions (large Pe) and non-spherical particles such as fibres.The effect of asymmetric flow after a bifurcation at inlet is also of great interest.Particle deposition due to cartilaginous rings in the trachea will also be considered.The flow in the trachea is in general turbulent and the present results which are valid for laminar flow can therefore then only be applied for very weak breathing conditions.Here it is of interest to investigate whether the separation due to the rings that occur for the laminar case will still be present in the turbulent case, since it is a well known phenomenon in fluid mechanics that turbulence resists the onset of separation.
In recent experiments (Åkerstedt et al. [15]) on carbon nanotubes it has been seen that the particles are electrically charged, therefore studies of deposition due to charged particles will also be considered.

Figure 1 .
Figure 1.Axisymmetric geometry of a simple cartilaginous ring configuration.

Figure 3 .
Figure 3. Local deposition rates for finite and zero value of .To connect the results with the ring structure, the shape of the rings f(x) are also shown.

Figure 4 .
Figure 4. Cumulative deposition P along the wall as a function of position X as obtained from (50) for different ringamplitudes.
minimum g(x) and for large x we have zero wall shear stress

Figure 5 .
Figure 5.Comparison between numerical and analytical solutions for smooth pipe.

Figure 6 .
Figure 6.Comparison with numerical and analytical results for a ring amplitude = 0.05.

Figure 7 .
Figure 7. Cumulative deposition P along the wall as a function.

Figure 8
Figure 8 Concentration distribution and streamlines for upper)and (lower).

Figure 10 .
Figure 10.Comparison of different analytical results (55-56) with numerical solution for a smooth pipe with developing velocity profile.

Figure 11 .
Figure 11.Cumulative deposition for different ring amplitudes for the case of developing flow.

Figure 12
Figure 12Cartilaginous ring configuration of Martonen type.

Figure 13 .
Figure 13.Cumulative deposition for the ring configuration by Martonen et al [12].