## Abstract

The paper discusses two *Full Lagrangian* methods for calculating the particle concentration and velocity fields in dilute gas-particle flows. The methods are mathematically similar but crucially different in application. By examining the analytical solution for two-phase stagnation-point flow, it is shown that Osiptsov's method is more general than that of Fernandez de la Mora and Rosner. In Osiptsov's method, the Jacobian of the Eulerian–Lagrangian transformation is computed by integration along particle pathlines. The particle concentration is then obtained algebraically from the Lagrangian form of the particle continuity equation. It is shown that the correct specification of the initial conditions is non-trivial and of vital importance. A technique to alleviate problems of mathematical ‘stiffness’ at small Stokes numbers is also described. Full Lagrangian methods require knowledge of the fluid velocity gradient field and, if the carrier flowfield is calculated numerically, differentiation of a ‘noisy’ field can result in serious errors. The paper describes a method for reducing these errors. The incompressible, inviscid flow over a cylinder provides a useful test case for validation and the Osiptsov method proves its worth by revealing a region, hitherto unknown, of crossing particle pathlines in the mathematical solution. Crossing pathlines and their relationship to Robinson's integral are then discussed, and calculations of particle flow through a turbine cascade at high Mach numbers are presented to illustrate the engineering potential of the method.

## 1. Introduction

In this paper, the term *dilute gas-particle flow* refers to those lightly laden particle flows where the particle motion is controlled by that of the carrier gas but not *vice versa* (one-way coupling), and where collisions between particles are unimportant. Examples of such aerosol-type flows abound in mechanical engineering (transport of solid combustion products, fouling of heat transfer equipment, particle and droplet deposition on turbine blades), chemical engineering (transport of chemically created aerosols), environmental science (transport and fallout of airborne pollutants), medical engineering (drug delivery to the lungs), agricultural engineering (crop spraying) and so on. Indeed, the claim that lightly laden particle flows represent the norm rather than the exception is hardly an exaggeration.

If the gas-particle flow is dilute in the sense described above, the gas flowfield is not affected by the particles and can be established by an independent calculation. This paper does not address the calculation of the gas flowfield, and the assumption is made that a suitable flow solution already exists. For general engineering calculations this would probably involve running a Reynolds-averaged Navier–Stokes, computational fluid dynamics (RANS CFD) package with built-in turbulence models. The time-mean particle velocity field can then be obtained using a Lagrangian method for integrating the particle equations of motion along the time-mean particle pathlines (often called *trajectories*), subject to the known gas velocity field. Lagrangian methods are easy to code and are particularly suitable for those cases where the particle inertia is high so that the particles are unresponsive to the fluid turbulence. When particle turbulence is important, the Lagrangian approach can still be used in conjunction with a more detailed gas turbulence model in ‘random-walk mode’, as described, for example, by Kallio & Reeks (1989). Lagrangian methods are also invariably used for gas-particle direct numerical simulations (DNS), where the objective is to study the particle behaviour in a numerically generated turbulent flow with all length-scales resolved (Squires & Eaton 1991). This is currently an important and popular area of research as it provides information about particle behaviour that would be very hard, if not impossible, to obtain by experimental measurement (for example, see the work on mixing and demixing of particles by turbulent structures by Wang & Maxey (1993)).

In the Lagrangian approach, the particle equation of motion is integrated along discrete particle pathlines (Crowe *et al*. 1996). The particle velocity field is then obtained by repeated application of this procedure for a large number of pathlines, which together encompass the region of interest. For high inertia particles that do not respond significantly to the gas-phase turbulence, surprisingly few pathline calculations are required to establish with good accuracy the particle velocity fields in quite complex three-dimensional flows by interpolation from the pathline solutions to a fixed Eulerian grid. Such calculations require only modest computational expense and are now commonplace.

In some applications, knowledge of the particle velocity field suffices, but in others it is necessary to know how the particle concentration varies along pathlines and also through the flowfield in general. Examples include the prediction of particle deposition rates on solid surfaces (Marchioli *et al*. 2003) and also DNS calculations to investigate the preferential accumulation of particles in turbulent structures (Marchioli & Soldati 2002). Both these areas are currently receiving considerable research effort aimed at providing fundamental understanding of the behaviour of small particles in complex flowfields.

Unlike the particle velocity field, however, accurate calculation of the particle concentration field is notoriously difficult. The traditional method (Crowe 1982) involves averaging the particle residence time for all pathlines crossing a small computational cell surrounding the point of interest. Unfortunately, extremely large numbers of pathlines are required to cover the region of interest and to reduce the errors to acceptable levels. For example, to reduce the local error to less than 0.1% (a very modest target in CFD), *each* cell in a two-dimensional calculation must be intersected by more than 10 000 pathlines. Computational costs are excessive and three-dimensional calculations are more-or-less precluded. Aside from anything else, the method lacks elegance.

Although not well known, there exists a superior approach, based on integrating equations for the particle concentration along individual pathlines. The origins of this *Full Lagrangian* method are obscure, but the first overt derivation of a relevant set of equations was given in an appendix to a paper by Fernandez de la Mora & Rosner (1981). These authors did not exploit the method further, but a decade later, Geller *et al*. (1993) used the same equations to compute the particle concentration fields in some inviscid, incompressible carrier flows that could be represented analytically. In Russia, a related line of development is associated with the names of Osiptsov (1998) and Tarasova & Tsirkunov (2000; see Tsirkunov 2001).

An examination of the Full Lagrangian approach suggests that it has the potential, quite literally, to revolutionize gas-particle flow calculations when a knowledge of the particle concentration field is required. This is because (i) it holds the promise of reducing CPU requirements by orders of magnitude (thus allowing the calculation of much more complex three-dimensional flows), (ii) it is likely to be more accurate (particularly in resolving the fine details of intricate flowfields), (iii) it can deal with steady and unsteady flows without change of algorithm and (iv) the mathematical formulation, aside from its elegance, provides a sound framework for interpreting the physics of the particle behaviour. Nevertheless, despite these very substantial potential advantages, the method has barely been exploited and there are many aspects of its formulation and practical application that require careful investigation. That is the purpose of this paper.

The paper starts by comparing the methods of Fernandez de la Mora & Rosner (1981) and Osiptsov (1998). This is done by considering a two-dimensional gas-particle stagnation-point flow with crossing particle pathlines for which an analytical solution exists. Superficially, the two methods appear rather similar, but close examination of the mathematics highlights a crucially important difference. Fernandez de la Mora and Rosner (FMR) viewed the particle mass conservation equation essentially as an Eulerian equation which could, with advantage, be solved along Lagrangian pathlines by working with the substantive derivative. Osiptsov's method, on the other hand, is firmly rooted in Lagrangian theory, the particle concentration being obtained from the Lagrangian form of the mass conservation equation by computing the change in volume of an element of ‘particle fluid’ along its pathline. This, it emerges, allows the Osiptsov, but not the FMR, method to handle certain types of singularity where the particle concentration becomes infinite (in a mathematical sense). This is a matter of pivotal importance for both CFD and DNS applications.

The paper then deals with three special considerations. The first is a semi-analytical integration method for reducing the computational time for very small particles. This is important because, as the particle inertia tends to zero, so too does the gas-particle velocity slip, and the governing equations become mathematically stiff. Secondly, experience with the analytical stagnation-point solution highlights the vital importance of setting the correct initial conditions for the particle concentration calculation. This is a non-trivial problem and is discussed fully in the text. Finally, a feature of any Full Lagrangian method is the need to know the velocity *gradient* field of the carrier fluid. Osiptsov (1998) and Geller *et al*. (1993) only applied their methods to carrier flowfields that could be described analytically, and which therefore had smooth continuous gradients. As yet, the Full Lagrangian method has not been exploited with numerically generated flowfields, which are inherently ‘noisy’. Unless care is taken, numerical differentiation to obtain velocity gradients can result in excessive cumulative errors in particle concentration. A procedure to avoid this problem is described.

The paper then turns to the validation and practical implementation of the Osiptsov method. The incompressible, inviscid flow over a cylinder is a useful test case because the velocity gradient field can be obtained analytically as well as by CFD techniques. Particle calculations using the Osiptsov approach demonstrate that the use of numerically generated flowfields do not necessarily lead to significant loss of accuracy. The method also proves its worth by revealing a small region where particle pathlines cross and give rise to infinite particle concentrations (in the mathematical sense). This hitherto unnoted phenomenon in cylinder-flow leads to a consideration of crossing particle pathlines and their relationship to an intriguing analytical result of Robinson (1956). Finally, calculations of the particle concentration field in a compressible gas flow through a turbine cascade illustrate the potential of the method for use in more complex engineering situations.

## 2. Two ways of computing the particle concentration

The discussion is restricted to dilute gas-particle flows with spherical mono-dispersed particles satisfying the usual ‘dusty-gas’ approximations. Thus, it is assumed that the particles are driven by the carrier gas but not vice versa, and that collisions between particles can be neglected. For simplicity, it is assumed that the only significant force acting on the particles is the steady-state viscous drag. Other forces such as those due to gravity, electrophoresis and thermophoresis can easily be added to the momentum equations if they are important.

### (a) The method of Fernandez de la Mora and Rosner

For the conditions specified above, the conservation equations of mass and momentum for the ‘particle fluid’ can be written in Eulerian form for plane two-dimensional flow as(2.1a)(2.1b)(2.1c)where *ρ*_{p} is the particle mass density (the product of concentration and particle mass), and (*V*_{x}, *V*_{y}) and (*U*_{x}, *U*_{y}) are the *x*- and *y*-components of the particle and gas velocity, respectively. In this Eulerian representation, *ρ*_{p}, *V*_{x}, *V*_{y}, *U*_{x} and *U*_{y} are all considered to be functions of the space and time coordinates *x*, *y* and *t*. The differential operator D/D*t*_{p}=∂/∂*t*+*V*_{x}∂/∂*x*+*V*_{y}∂/∂*y* is the substantive derivative following an element of ‘particle fluid’, which (because diffusion is absent) always comprises the same particles. *β* is the reciprocal of the particle relaxation time, which, for spherical particles of diameter *d*_{p} and material density *ρ*_{mat} subjected to Stokes drag by a gas of dynamic viscosity *μ*_{g} would be given by(2.2)For simplicity, *β* is assumed to be constant throughout the flowfield.

Equations (2.1*b*) and (2.1*c*) can easily be integrated numerically along particle pathlines, but equation (2.1*a*) proves more problematic because of the difficulty in evaluating the divergence of the particle velocity. The solution offered by Fernandez de la Mora & Rosner (1981) was to derive additional equations for the substantive derivatives of ∂*V*_{x}/∂*x*, ∂*V*_{x}/∂*y*, ∂*V*_{y}/∂*x* and ∂*V*_{y}/∂*y* by differentiating equations (2.1*b*) and (2.1*c*) first by *x* at constant *y* and then by *y* at constant *x*. The extra four equations generated are(2.3a)(2.3b)(2.3c)(2.3d)Equations (2.1*a–c*) and (2.3*a–d*) can be integrated along particle pathlines, thus allowing the calculation of *ρ*_{p} alongside *V*_{x} and *V*_{y}. Initial conditions are required for *ρ*_{p}, *V*_{x}, *V*_{y}, ∂*V*_{x}/∂*x*, ∂*V*_{y}/∂*y*, ∂*V*_{y}/∂*x* and ∂*V*_{y}/∂*y*. Values of the fluid velocity components *U*_{x} and *U*_{y}, and the fluid velocity gradients ∂*U*_{x}/∂*x*, ∂*U*_{x}/∂*y*, ∂*U*_{y}/∂*x* and ∂*U*_{y}/∂*y* are also required at each integration point on a pathline. Assuming these are available, numerical integration is straightforward except at points where the divergence of the particle velocity field is undefined because ∂*V*_{x}/∂*x* or ∂*V*_{v}/∂*y*, or both, become infinite. As explained below, however, this can provide a very undesirable constraint on the usefulness of the method.

### (b) The method of Osiptsov

In Lagrangian form, the conservation equations of mass and momentum for an element of particle fluid are written(2.4a)(2.4b)(2.4c)In this *Lagrangian* representation, *ρ*_{p}, *V*_{x}, *V*_{y}, *U*_{x} and *U*_{y} are all considered to be functions of the Lagrangian independent variables *x*_{p0}, *y*_{p0} and *τ*. Here, *τ* is the time and the space variables *x*_{p0} and *y*_{p0} pick out a particular pathline by the fact that the particle element passes through the point (*x*_{p0}, *y*_{p0}) at *τ*=0. The partial derivatives ∂*V*_{x}/∂*τ* and ∂*V*_{y}/∂*τ* are evaluated at constant *x*_{p0} and *y*_{p0}, and have the same numerical values as the substantive derivatives D*V*_{x}/D*t*_{p} and D*V*_{y}/D*t*_{p} in equations (2.1*b*) and (2.1*c*). *U*_{x} and *U*_{y} in equations (2.4*b*) and (2.4*c*) are the components of the gas velocity ‘seen’ by the particles at time *τ*.

In equation (2.4*a*), *ρ*_{p0} is the mass density of the particle fluid element at *τ*=0. At this time the element is rectangular of area *δx*_{p0}*δy*_{p0} (i.e. two-dimensional volume) as shown in figure 1. At a later time *τ*, it has become lozenge-shaped and has area |*J*|*δx*_{p0}*δy*_{p0}. Just as in single-phase Lagrangian fluid mechanics, it can be shown that, in the limit as *δx*_{p0} and *δy*_{p0} both tend to zero,(2.5a)where *J*_{xx}, *J*_{xy}, *J*_{yx} and *J*_{yy} are given by(2.5b)In equations (2.5*b*), the precise meaning of each partial derivative has been made clear by inclusion of subscripts which emphasize exactly which variables are to be held constant. The geometrical interpretation of these derivatives is shown in figure 1. *J* can also be interpreted as the Jacobian of the Eulerian–Lagrangian transformation. If particle pathlines cross, the Jacobian changes sign, although the area of the element remains positive. That is the reason for requiring the absolute value of *J* in equation (2.4*a*).

Four new variables *w*_{xx}, *w*_{xy}, *w*_{yx} and *w*_{yy} are now defined by the relationships(2.6a)(2.6b)(2.6c)(2.6d)Differentiating equations (2.6*a–d*) with respect to *τ* and introducing equations (2.4*b*) and (2.4*c*) then gives(2.7a)(2.7b)(2.7c)(2.7d)The fluid velocity components *U*_{x} and *U*_{y} are normally only known in Eulerian coordinates, but the required Lagrangian derivatives can be obtained using the chain rule. Thus,(2.8a)

(2.8b)

(2.8c)

(2.8d)

Equations (2.4*a–c*), (2.6*a–d*) and (2.8*a–d*) can be integrated numerically along pathlines to give *ρ*_{p}, *V*_{x} and *V*_{y}. Initial conditions are required for *ρ*_{p}, *V*_{x}, *V*_{y}, *J*_{xx}, *J*_{xy}, *J*_{yx}, *J*_{yy}, *w*_{xx}, *w*_{xy}, *w*_{yx} and *w*_{yy}. Also required are values of *U*_{x}, *U*_{y}, ∂*U*_{x}/∂*x*, ∂*U*_{x}/∂*y*, ∂*U*_{y}/∂*x* and ∂*U*_{y}/∂*y* at each integration point on the pathline. Although the Osiptsov method involves more differential equations than that of FMR, the computational requirements are very similar. Crucially, however, in Osiptsov's method, *J* is treated as the dependent variable and is found by integrating the *differential* equations (2.6*a–d*) and (2.8*a–d*). The particle concentration *ρ*_{p} is then obtained by a ‘post-processing’ application of equation (2.4*a*), which is *algebraic*. In the FMR method, the dependent variable is *ρ*_{p}, which is obtained directly by integrating the *differential* equation (2.1*a*). This gives rise to numerical difficulties whenever *ρ*_{p} becomes infinite. On the other hand, Osiptsov's method can handle singularities in *ρ*_{p}, because at these points, *J* passes smoothly through zero.

An interesting discussion of different types of mathematical singularity in gas-particle flow is given by Osiptsov (1984). It should, of course, be appreciated that in real flows infinite values of *ρ*_{p} can never occur because the effects of finite particle size and inter-particle collisions become important at high particle concentrations. Nevertheless, the equations for dilute gas-particle flow remain an important modelling approximation (just as the Euler equations with singularities in the form of vortex lines are an important approximate form of the Navier–Stokes equations in single-phase flow), and provide much physical insight as well as being the basis for an indispensable computational technique.

### (c) Particle flow near a stagnation-point

The flow of particles in an incompressible, inviscid carrier gas near a stagnation-point provides one of the few known analytical solutions of the particle equations of motion. It is particularly useful for assessing different calculation schemes because, depending on the Stokes number, quite different types of particle behaviour are exhibited. This flowfield is now used to illustrate the differences in the two methods for calculating the particle concentration.

The gas flowfield is given by the velocity potential *ϕ*_{g}=*A*(*y*^{2}−*x*^{2})/2 (with *A* a positive constant). The velocity components and their derivatives are then(2.9a)(2.9b)and the gas streamlines are shown in figure 2*a*. The gas flows enter the domain through the left and right boundaries and are turned so that they exit through the top and bottom boundaries. The stagnation-point is at the origin of coordinates in the centre (*x*=0, *y*=0).

The particle pathlines are obtained by substituting the relationships *V*_{x}=∂*x*_{p}/∂*τ* and *V*_{y}=∂*y*_{p}/∂*τ* together with equations (2.9*a*,*b*) into equations (2.4*b*) and (2.4*c*), thus giving two second-order differential equations for *x*_{p} and *y*_{p} in terms of the time *τ*:(2.10a)(2.10b)*St*=*A*/2*β* is the Stokes number, the ratio of the particle relaxation time 1/*β* to a flow time-scale 2/*A*.

Equations (2.10*a*,*b*) admit two different types of solution. For *St*<0.125 (subcritical flow), particle pathlines, starting far upstream of the stagnation-point with zero velocity slip, do not cross other particle pathlines nor the *y*-axis. This behaviour has been discussed fully by Fernandez de la Mora & Rosner (1981) and is illustrated in the lower diagram of figure 2*b*, which shows, for *St*=0.05, gas streamlines and particle pathlines (in the upper left quadrant only). The particle density *ρ*_{p} can be shown to be independent of *y*, and its variation with *x* is plotted in the upper diagram of figure 2*b*. Although not obvious from the graph, *ρ*_{p} becomes infinite on the *y*-axis (*x*=0).

In contrast, for *St*>0.125 (supercritical flow) particles oscillate about the *y*-axis with decaying amplitude and the pathlines cross. This behaviour can be seen in the lower diagram of figure 2*c*, which shows, for *St*=0.5, particles entering the upper left quadrant from the left, crossing into the upper right quadrant, and then oscillating around the *y*-axis with decreasing amplitude before exiting the domain through the top boundary. The variation of *ρ*_{p} along pathlines is more complicated than for subcritical flow. As can be seen by a careful study of the upper diagram of figure 2*c*, *ρ*_{p} becomes infinite at each ‘turning point’ where the *x*-component of particle velocity *V*_{x} is zero.

Instead of relying on the analytical solution, suppose instead that *V*_{x}, *V*_{y} and *ρ*_{p} are computed using a numerical stepwise integration along particle pathlines (for any specified initial conditions). For subcritical flows, either of the two Lagrangian methods discussed above could be used successfully. The only slight difficulty would occur on the stagnation pathline, where particles require an infinite time to reach the stagnation-point where *ρ*_{p}→∞.

The same is not true for supercritical flow, however. The analytical solution for the coordinates of the particle pathlines in supercritical flow, obtained from equations (2.10*a*,*b*) for arbitrary initial conditions (*x*_{p}=*x*_{p0}, *y*_{p}=*y*_{p0}, *V*_{x}=*V*_{x0}, *V*_{y}=*V*_{y0} at *τ*=0), is(2.11a)(2.11b)In equations (2.11*a*,*b*), (*x*_{p}, *y*_{p}) are the coordinates of a point on the particle pathline at time *τ*, *T*=*βτ*/2 is a non-dimensional time, *γ*=(8*St*+1)^{1/2} and *α*=(8*St*−1)^{1/2}.

For illustration, consider the particular solution corresponding to the initial conditions *V*_{x0}=*U*_{x0}=−*Ax*_{p0} (zero initial *x*-velocity slip) and *V*_{y0}=0 (zero initial particle *y*-velocity). From the discussion above, it is expected that the pathlines cross the *y-*axis and then oscillate around it with decaying amplitude (as in figure 2*c*). This behaviour is shown in a different way in figure 3*a*, where *x*_{p}/*x*_{p0} is plotted against non-dimensional time *T* (the solid line). The particle velocity components (*V*_{x}, *V*_{y}) for these initial conditions are given analytically by(2.12a)(2.12b)and the variation of *V*_{x}/*U*_{x0} with *T* is also shown in figure 3*a* (the dotted line). The condition *V*_{x}=0 corresponds to the ‘turning points’ of the pathlines and occurs when tan(*αT*)=−*α*.

Now consider the variation of the Jacobian with *T*. Because *V*_{x} is a function only of *x*, the component *J*_{xy} is zero everywhere and hence *J*=*J*_{xx} *J*_{yy}. Combining equations (2.6*a*) and (2.8*a*), and (2.6*d*) and (2.8*d*) with ∂*U*_{x}/∂*y*=∂*U*_{y}/∂*x*=0 shows that *J*_{xx} and *J*_{yy} satisfy the differential equations(2.13a)(2.13b)When solving equations (2.13*a*,*b*), it is most important to specify the initial conditions correctly. The conditions (*J*_{xx})_{0}=(*J*_{yy})_{0}=1 pose no difficulty, but the specification of (∂*J*_{xx}/∂*τ*)_{0} and (∂*J*_{yy}/∂*τ*)_{0} requires care. A full discussion is given in §3*b* below, where it is shown that(2.14a)(2.14b)The solution of equations (2.13*a*,*b*) with these initial conditions is then(2.15a)(2.15b)Equations (2.15*a*,*b*) can also be obtained by applying the definitions *J*_{xx}=∂*x*_{p}/∂*x*_{p0} and *J*_{yy}=∂*y*_{p}/∂*y*_{p0} to equations (2.11*a*,*b*) directly and then specializing to the particular solution. As shown in figure 3*a*, *J*_{xx} varies with *T* in the same way as *V*_{x}/*U*_{x0}.

Equations (2.15*a*,*b*) show that *J*_{yy}>0 for all *T*, but *J*_{xx}=0 when tan(*αT*)=−*α*. Bearing in mind that *J*=*J*_{xx}*J*_{yy}, it can be seen from equation (2.4*a*) and figure 3*b* that *ρ*_{p} becomes infinite when tan(*αT*)=−*α*. This, of course, corresponds to the condition *V*_{x}=0 on the particle pathlines. Equations (2.15*a*,*b*) show that *J*_{xx} and *J*_{yy} are well-behaved at all times. A stepwise numerical integration along trajectories using the Osiptsov method would have no difficulty in the vicinity of the turning points defined by tan(*αT*)=−*α*, *even though* *ρ*_{p} *becomes unboundedly large*. The success of the Osiptsov procedure has, indeed, been confirmed by numerical integration of the equations for various supercritical Stokes numbers. The results have not been plotted in figure 3 because the curves are indistinguishable from the analytical solutions.

Turning now to the FMR method, equations (2.3*a*) and (2.3*d*) for stagnation-point flow become(2.16a)(2.16b)With the initial conditions ∂*V*_{x}/∂*x*=∂*V*_{y}/∂*y*=0, the solution of equations (2.16*a*,*b*) is(2.17a)(2.17b)Equations (2.17*a*,*b*) can be substituted into equation (2.1*a*) and integrated to give an expression for *ρ*_{p} that is consistent with that for *J* obtained from equations (2.15*a*,*b*). Unfortunately, however, problems occur if a numerical solution is attempted. Although ∂*V*_{y}/∂*y* is well behaved, equation (2.17*a*) shows that ∂*V*_{x}/∂*x* becomes infinite (along with *ρ*_{p}) when tan(*αT*)=−*α*. This is also shown by the curve for the divergence of the particle velocity (∂*V*_{x}/∂*x*+∂*V*_{y}/∂*y*) in figure 3*b*.

This instructive example is extremely useful for obtaining a ‘feel’ for the way the equations behave. It shows that, by solving for the Jacobian, the Osiptsov method has the advantage over the FMR approach of being able to handle those particle concentration singularities where the Jacobian changes sign smoothly. This is crucially important for a numerical scheme, particularly for DNS gas-particle calculations where the nature of the carrier flowfield is such that mathematical singularities in particle concentration can occur randomly. The FMR method would simply not be able to cope with such situations. Hence, because of its wider range of applicability, Osiptsov's method has been used for all subsequent calculations presented in this paper.

## 3. Special considerations for the Osiptsov method

### (a) Treatment for small Stokes numbers

Equations (2.4*a–c*), (2.6*a–d*) and (2.8*a–d*) of the Osiptsov method can be integrated numerically along particle pathlines using any convenient algorithm for simultaneous ODEs. The authors have found that a simple predictor–corrector scheme usually gives good accuracy. If an explicit integration method is used, the integration time-step must be less than the particle relaxation time *β*^{−1} in order to ensure numerical stability. This restriction does not pose a problem at high Stokes numbers because the time-step is usually limited, for reasons of accuracy, by the requirement that the carrier flow conditions do not change significantly over a single integration step. At small Stokes numbers, however, the equations become mathematically ‘stiff’, very small time-steps are required and computational times become excessive. One method of overcoming this constraint is to use an approximate analytical integration of the equations over a time-interval long compared with *β*^{−1}, but short compared with the time-scale of the carrier flow changes. This approach is simpler and more transparent than traditional backward-difference numerical methods for stiff ODEs.

Consider first the particle equation of motion (2.4*b*) written in the form(3.1)Multiplying by e^{βτ} and integrating from *τ*=0 to *τ*=Δ*τ* assuming constant ∂*U*_{x}/∂*τ* gives(3.2)where the overbar on ∂*U*_{x}/∂*τ* denotes an average value over the time-step. The first term on the right-hand side represents the exponential decay with time-constant *β*^{−1} of the residual velocity slip existing at *τ*=0. The second term represents the exponential approach to the local ‘steady-state’ slip velocity . Equation (3.2) can be implemented using a similar algorithm as for equation (2.4*b*), but now there is no stability restriction and Δ*τ* is only limited by the accuracy of the assumption that ∂*U*_{x}/∂*τ* remains constant. Equation (2.4*c*) can be handled in exactly the same way.

A similar approach can be used to integrate out the stability restriction in equations (2.8*a–d*). Multiplying equation (2.8*a*) by e^{βτ} and integrating analytically assuming constant *J*_{xx}∂*U*_{x}/∂*x* and *J*_{yx}∂*U*_{x}/∂*y* gives(3.3)A further integration gives an expression for *J*_{xx} to replace equation (2.6*a*):(3.4)Equations (2.8*b–d*) can be handled in similar ways.

### (b) Initial conditions for the Jacobian calculation

In his 1998 paper, Osiptsov did not discuss the specification of the initial conditions for the Jacobian calculation, but this problem is non-trivial and requires careful consideration.

Figure 4 shows an element of ‘particle fluid’ at *τ*=0. This element is rectangular irrespective of the local pathline directions. Hence, from the definitions (2.5*b*),(3.5a)(3.5b)The derivatives (*w*_{xx})_{0}, (*w*_{xy})_{0}, (*w*_{yx})_{0} and (*w*_{yy})_{0} are more difficult to specify unless the calculation starts in a region where the carrier flow is uniform with streamlines parallel to the *x*-axis, and the particle velocity slip is zero. In such cases, it follows from equations (2.6*a–d*) that(3.6)

In more general cases, the carrier flow may be non-uniform and the initial velocity slip non-zero. A common example (see figure 4) is when the particles are released in a non-uniform (but well-defined) carrier velocity field with the particle velocity variation at *τ*=0 specified along a particular line *x*_{p0}=*x*_{0}=*constant* in the form *V*_{x0}=*V*_{x0}(*y*_{p0}) and *V*_{y0}=*V*_{y0}(*y*_{p0}). At *τ*=0, the particle velocities at points 1 and 3 in figure 4 are defined directly by the initial conditions, and this is sufficient to specify (*w*_{xy})_{0} and (*w*_{yy})_{0}. Thus, from equations (2.6*b*) and (2.6*d*),(3.7a)(3.7b)At time *τ*=0, the particle velocity at point 2 in figure 4 is fixed, not just by the particle velocities specified along *x*_{p0}=*x*_{0}, but also by the particle velocity variation along the pathline 1′→2. Starting from the definitions, equations (2.6*a*) and (2.6*c*), some consideration shows that the correct expressions for the remaining components (*w*_{xx})_{0} and (*w*_{yx})_{0} are(3.8a)(3.8b)The first term on each right-hand side represents the change from 1→1′ and the second term the change from 1′→2.

As an example, consider the stagnation-point flow discussed previously. The initial conditions *V*_{x0}=*U*_{x0} and *V*_{y0}=0 lead straightforwardly to (*w*_{xy})_{0}=0 [from equation (3.7a) because ∂*U*_{x}/∂*y*=0 everywhere] and (*w*_{yy})_{0}=0 [from equation (3.7*b*) because *V*_{y0}=0]. Also, (*w*_{xx})_{0}=0 [from equation (3.8*a*) because *V*_{y0}=0 and ], and (*w*_{yx})_{0}=*βU*_{y0}/*U*_{x0} [from equation (3.8*b*) because *V*_{y0}=0 and hence (∂*V*_{y}/∂*τ*)_{0}=*βU*_{y0}].

It should be noted that (*w*_{yx})_{0} (which is non-zero) was not required for the analysis of §2*c* because the special nature of stagnation-point flow was acknowledged *a priori* by the assumption (which turns out to be correct) that *V*_{x} and *V*_{y} are functions of *y* only and hence that *J*_{xy}=0. Had this assumption not been made, and had it been assumed erroneously that (*w*_{yx})_{0}=0, an incorrect expression for *J* would have resulted in which the singularities in *ρ*_{p} occur, not at the turning points where *V*_{x}=0, but on the *y*-axis. This example illustrates the importance of setting the initial conditions correctly and how errors in the Jacobian or its derivatives, unlike errors in particle velocity, do not decay exponentially with time constant 1/*β*. Such issues only really become apparent, however, if the mathematics for the stagnation flow are worked through first with the correct and then with incorrect starting assumptions.

### (c) Calculation of the fluid velocity gradients

Integration of equations (2.8*a–d*) requires a knowledge of the local fluid velocity gradients ∂*U*_{x}/∂*x*, ∂*U*_{x}/∂*y*, ∂*U*_{y}/∂*x* and ∂*U*_{y}/∂*y*. If the carrier flowfield is represented analytically (as in the work of Osiptsov (1984) and Geller *et al*. (1993)), these can be obtained by analytical differentiation and will automatically vary smoothly. If, however, the Full Lagrangian method is to prove of real practical value, then it must operate accurately with numerically generated flowfields. Currently, one of the most popular types of CFD solver is the finite-volume time-marching method, and it is known that particle *velocity* fields can be computed accurately, based on solutions of this type. What is not known is whether or not the numerical differentiation of a ‘noisy’ CFD solution will generate unacceptable errors for Lagrangian particle concentration calculations given that (as shown in §3*b*) errors incurred along pathlines do not decay and are therefore cumulative.

Numerical differentiation enhances errors so, to provide as smooth a velocity gradient field as possible, the derivatives were obtained using a technique based on numerical integration. This method, which is sometimes used to obtain derivatives in single-phase Navier–Stokes solvers, can be applied to any arbitrary non-uniform grid. It is based on the ‘gradient theorem’ (a special case of Gauss's theorem) and is expressed by(3.9)where *vol* is a finite volume enclosed by a surface *S* and *ϕ* is any scalar function of position. Application of equation (3.9) in discretized form to the two-dimensional quadrilateral cell shown in figure 5 gives(3.10)On the left-hand side, *L*_{x} and *L*_{y} are the projections of a side normal to the *x*- and *y*-directions and *ϕ* represents a mean value along the side. On the right-hand side, ∂*ϕ*/∂*x* and ∂*ϕ*/∂*y* are mean values of the derivatives, averaged over the cell of area *A*. Setting *ϕ*=*U*_{x} and then *ϕ*=*U*_{y} and separating into components gives the four scalar equations(3.11a)(3.11b)Equations (3.11*a*,*b*) can be used to calculate the four required velocity gradients from the CFD-generated velocity field. Values along particle pathlines are then obtained by interpolation. The method can easily be adapted for either cell-centred or cell-cornered storage schemes.

## 4. Gas-particle flow over a circular cylinder

### (a) Calculation of the gas flowfield

Inviscid, incompressible two-dimensional flow over a cylinder is unphysical because there are no boundary layers, but it is still a useful test case because of the existence of an analytical solution for the velocity field. The velocity potential *ϕ*_{g} is given by(4.1)where *U*_{∞} is the oncoming flow velocity and *R* is the cylinder radius. Analytical expressions for the velocity components and their derivatives can be obtained by differentiating equation (4.1).

The velocity field was also obtained from an inviscid CFD time-marching code written for ductflows. In order to model the flow over an isolated cylinder, the duct height was made very large in comparison to the cylinder diameter. The near-cylinder CFD grid is shown in figure 6, but the actual extent in the *y*-direction was five times greater. Figure 7 shows a comparison between the analytical and CFD-based calculations of ∂*U*_{x}/∂*x* and ∂*U*_{x}/∂*y*. Similar agreement was obtained for ∂*U*_{y}/∂*x* and ∂*U*_{y}/∂*y*. Although the main features are captured well, it can be seen that the agreement in some regions is not as good as might be achieved with a finer computational mesh. Nevertheless, the CFD solution was deliberately not refined because the objective was to establish if the particle concentration could be computed with acceptable accuracy using a flowfield generated as in a typical practical application.

### (b) Calculation of the particle concentration field

Particles were computationally injected along the line *x*_{p0}=−3*R* with uniform mass density *ρ*_{p0}, zero *x*-direction velocity slip *V*_{x0}=*U*_{x0}, and zero *y*-direction velocity *V*_{y0}=0. The initial conditions for equations (2.6*a–d*) and (2.8*a–d*) are given by equations (3.5*a*,*b*), (3.7*a*,*b*) and (3.8*a*,*b*) as(4.2a)(4.2b)The Osiptsov equations were then integrated numerically, the fluid velocity gradients being obtained either analytically using equation (4.1) or from the CFD solution using equations (3.11*a*,*b*).

The results are presented in the form of contours of normalized particle concentration *ρ*_{p}/*ρ*_{p0}. Contouring directly from Lagrangian data is not to be recommended so the results were first interpolated onto the CFD mesh of quasi-streamlines and quasi-orthogonals using the method of Geller *et al*. (1993). In brief, an initial interpolation along pathlines gave a set of off-mesh values on each quasi-orthogonal, and a second interpolation along the quasi-orthogonals gave a set of on-mesh values, suitable for contouring.

As a further comparison, the concentration field was also found using the traditional ‘counting’ method of Crowe (1982). This involved computing a very large number of pathlines in order that every Eulerian cell was intersected by a ‘substantial’ number. For each cell, the sum of the transit times of all intersecting pathlines was recorded, the particle concentration then being proportional to this sum.

### (c) Results and conclusions

The Stokes number for the cylinder calculations is defined by *St*=*U*_{∞}/*Rβ*, this being numerically compatible with the definition *St*=*A*/2*β* used for the stagnation-point analysis. The steady-flow of a dust-laden incompressible, inviscid gas around a sphere at small Stokes numbers has been analysed by Michael (1968) and the flow around a cylinder has many features in common. Thus, when *St* is less than a critical value *St*_{crit}=0.125, particle pathlines do not intersect the cylinder surface. For *St*>*St*_{crit} particles reach the stagnation-point with finite velocity and deposition may occur.

Figure 8 shows the normalized particle concentration fields for *St*=10, 1, 0.1 and 0.01, obtained using the Osiptsov method. In each case, 500 pathlines were initiated and the integration time-step was chosen so that about five steps were needed for the particles to cross the smallest cell of figure 6. The contour plots on the left of figure 8 refer to calculations using equations (3.11*a*,*b*) and the CFD solution, and those on the right to calculations based on equation (4.1).

For *St*=10 (figure 8*a*) the particles have high inertia and their pathlines are almost straight. Particles launched directly upstream of the cylinder deposit on the surface (a perfectly absorbing boundary is assumed) and a particle-free zone forms behind the cylinder. Both calculations capture this zone accurately, but there are differences in the concentration fields elsewhere. Closer examination shows, however, that these variations are quite small. Figure 8*b–d* show the calculations for *St*=1, 0.1 and 0.01, respectively. In all cases, the CFD-based solution is very similar to that obtained using the analytical equations.

A more stringent test of accuracy is obtained by examining the variation of *ρ*_{p}/*ρ*_{p0} along the outflow boundary as shown in figure 9. The results of four different calculations are given on each graph. The solid and dotted lines correspond to ‘Osiptsov CFD’ and ‘Osiptsov analytical’ using 500 pathlines as in figure 8. The squares and circles correspond to ‘traditional CFD’ and ‘traditional analytical’ calculations. These calculations required 10 000 pathlines.

For *St*=10 (figure 9*a*), *ρ*_{p}/*ρ*_{p0} at outlet is close to unity and all four calculations perform well. For *St*=1 and 0.1 (figure 9*b*,*c*), the ‘Osiptsov CFD’ and ‘Osiptsov analytical’ results show similar maximum values at the edge of the particle-free zone, albeit at slightly different locations. For *St*=0.01 (figure 9*d*), both Osiptsov calculations are in close agreement.

Considering the ‘traditional’ calculations of figure 9*b–d*, it can be seen that the results based on the analytical flowfield are in closer agreement with the Osiptsov calculations than those based on the CFD solution. In order to locate the source of the discrepancies, further calculations were undertaken and the following findings emerged:

For the Osiptsov calculations, using more than 500 pathlines had no effect. Indeed, as few as 200 pathlines gave essentially the same results.

The errors incurred with a CFD-based Osiptsov method are typically as shown in figure 9. These might be reduced by using a more refined CFD calculation or, better still, a coupled CFD-pathline calculation with adaptive grid refinement in regions of high concentration gradient.

For the traditional calculations, using more than 10 000 pathlines with the same Eulerian grid had no effect. However, refining the grid size by a factor of 2 in the

*y*-direction and increasing the number of pathlines to 20 000 (to maintain the same number of intersections per cell) gave much closer agreement with the Osiptsov results.

In general, the traditional counting method is not to be recommended as it requires excessive pathline calculations and a fine Eulerian grid. A two-dimensional Osiptsov calculation requires the solution of 10 ODEs compared with 2 for the traditional method. Nevertheless, the former is still some 20 times more efficient than the latter. In three dimensions the improvement would be even more dramatic. Indeed, the traditional method would not really be viable.

## 5. Crossing particle pathlines

### (a) Relationship with Robinson's analysis

Crossing particle pathlines have already been discussed in connection with the stagnation-point flow. They also occur in inviscid gas-particle flow over a cylinder, although this does not seem to have been reported in the literature. The phenomenon bears a relationship to the intriguing analysis of Robinson (1956), which is now discussed.

Robinson's analysis concerns the behaviour of particles in carrier flows, which can be represented by a velocity potential *ϕ*_{g}. Both the stagnation-point and cylinder flows discussed above fall into this category. By considering the particle equation of motion, Robinson showed that Kelvin's ‘circulation theorem’ holds for the particle as well as the carrier flowfield. Thus, if the particles start from far upstream with zero velocity slip, the circulation around a closed loop based on either the fluid or the particle velocity is zero, and will remain zero (for both particles and fluid) as the loops are convected in their different ways through the flowfield. Given that Kelvin's theorem holds for the particles, Robinson argued that a velocity potential *ϕ*_{p} for the particle flowfield can be defined. The remarkable result then followed that the particle concentration cannot decrease along a pathline, D*ρ*_{p}/D*t*_{p}≥0.

Robinson's work is not widely known in the two-phase flow community, although it has been used to effect by Fernandez de la Mora & Rosner (1981) in their analysis of subcritical stagnation-point flow. Sometimes, however, it is invoked incorrectly. For example, the unqualified statement by Geller *et al*. (1993) that ‘the (particle) concentration cannot decrease along a particle trajectory in potential flow’ is misleading and incorrect.

The validity of Robinson's analysis is related to the existence of crossing particle pathlines. Thus, if Robinson's result were to hold for all potential carrier flows (subject to the far upstream requirement that *ϕ*_{p}→*ϕ*_{g}), then it would follow that the crossing of particle pathlines in any potential carrier flow was forbidden (because *ϕ*_{p} would be indeterminate at the crossing points). That pathlines in a potential carrier flowfield can indeed cross is easily seen from the equations for supercritical stagnation-point flow. Interestingly, FMR (using Robinson's approach) found that the velocity potential *ϕ*_{p} did not exist when the Stokes number exceeded the critical value.

The problem in Robinson's analysis can be traced to the assumption that Kelvin's circulation theorem automatically implies the existence of a velocity potential. Figure 10 shows how Kelvin's theorem continues to hold even if particle pathlines cross. Thus, if the circulation around loop (1234) is zero, it will remain zero, not only around loop (1′2′3′4′), but also around loop (1″2″3″4″). The latter case obviously precludes the possibility of defining a velocity potential and the subsequent analysis leading to D*ρ*_{p}/D*t*_{p}≥0 is no longer valid. The situation does not arise in fluid flows because fluid streamlines can never cross.

### (b) Infinite particle concentration in cylinder flow

Figure 11 shows an example of crossing particle pathlines in a cylinder flow. The Stokes number is 0.2, which is slightly supercritical. In order to aid interpretation, figure 11*a* is plotted with pathline 2 (used as a reference pathline) projected onto the abscissa. Particles are launched at *X*=*x*/*R*=0, the cylinder leading edge is at *X*=2 and the centre is at *X*=3.

Particle pathlines 1–8 are plotted so that the ordinate represents their displacement from pathline 2. Pathline 1 deviates downwards before impacting the cylinder. Pathlines 3–6 initially deviate upwards before ‘closing the gap’ and crossing the reference pathline and their neighbours one after the other. Pathlines 7 and 8 do not intersect their nearest neighbouring pathlines. Figure 11*b* has the same abscissa scale as figure 11*a*, but the ordinate scale is expanded by a factor of 20 to show the intersections in greater detail.

Figure 11*c* shows the variation of the Jacobian along each of the trajectories. Clearly, *J* passes smoothly through zero (and *ρ*_{p} becomes infinite) as each of pathlines 3–6 crosses its nearest neighbour. Subsequently, |*J*| increases and the particle concentration decreases (D*ρ*_{p}/D*t*_{p}<0) showing that Robinson's result is not applicable despite the potential nature of the carrier flowfield. Although Osiptsov's method has no difficulty, the FMR method would not be able to cope with this type of calculation because the divergence of the particle velocity field becomes infinite along with *ρ*_{p}. It is interesting to note that Geller *et al*. (1993) who used the FMR method to calculate a similar flow make no mention of this difficulty, or of the existence of lines of infinite concentration. A numerical investigation conducted by the present authors revealed that lines of infinite *ρ*_{p} exist throughout the range 0.2<*St*<2.0.

## 6. Gas-particle flow through a turbine cascade

A final, and more complex, example of the usefulness of the Osiptsov method is that of gas-particle flow through a turbine cascade. This choice of geometry reflects the authors' interest in particle transport in turbomachinery because of its importance in clean-coal technology for future generations of gas turbines. In a coal gasification plant, small coal ash particles that are not removed by the filter system enter the turbine and may deposit on the blades or cause erosion. In predicting the magnitude of these effects, it is important to be able to calculate both the particle concentration and velocity fields within the blade rows.

Unlike the stagnation-point and cylinder flows discussed previously, the cascade carrier flow was compressible with large variations in gas density (the exit Mach number was 1.12). It was computed using the well-known CFD flow solver developed by Denton (1992). This code solves the Reynolds-averaged Navier–Stokes equations with a mixing-length turbulence model and can handle mixed subsonic-supersonic flows. The particle density field was calculated using the Osiptsov method, and figure 12 shows the results for Stokes numbers *St* of 1.0 and 0.1. (In calculating *St*, the ‘flow time’ was based on the blade axial chord length and mean axial flow velocity.) As expected, a particle-free zone forms behind the blade in both cases, the low-inertia particles (*St*=0.1) responding more to the curvature of the streamlines.

It is difficult to assess the accuracy of concentration calculations in such a complex flow. However, for any compressible gas flow, a useful check can be made in the limit *St*→0 when particle pathlines and gas streamlines should coincide. The Eulerian continuity equations are(6.1)With zero velocity slip, D/D*t*_{p}=D/D*t*_{g}=D/D*t* and equation (6.1) can be integrated to give(6.2)This shows that, when suitably non-dimensionalized, the particle concentration field in the limit *St*→0 should be identical to the gas density field. A comparison between the two at a very low Stokes number therefore provides a sensitive test of computational accuracy. Figure 13 shows that the Osiptsov method performs well in this respect.

## 7. Conclusions

The paper has investigated the potential of two Full Lagrangian methods for calculating the particle concentration field in two-dimensional flows. It has been shown that the Osiptsov approach is preferable to that of Fernandez de la Mora & Rosner (1981) because it can handle certain types of particle concentration singularity. It is, however, of paramount importance that the initial conditions for the Jacobian and its derivatives be specified correctly. A special technique has been developed for flows at a low Stokes number that relieves the mathematical stiffness of the equations and allows solutions to be obtained without increase in CPU time. It has also been established that Osiptsov's method can be used with CFD-generated flowfields and a procedure for calculating the required fluid velocity derivatives from ‘noisy’ velocity fields without introducing unacceptable errors has been proposed. The possibility of intersecting particle pathlines has been discussed in the context of Robinson's analysis and it has been emphasized that, unlike fluid flows, the validity of Kelvin's circulation theorem for particles does not necessarily imply the existence of a velocity potential. It has also been shown that crossing trajectories and lines of infinite particle concentration occur in the mathematical solution of two-dimensional gas-particle cylinder flow. The fact that such phenomena have not been observed previously in what is a comparatively simple flowfield suggests that similar behaviour might be found in more complex flows.

Osiptov's method for calculating the particle concentration field has great potential for dramatic reductions in computational time and improvements in accuracy compared with the traditional approach. Particle concentration calculations in industrial three-dimensional flows, previously precluded by prohibitive computational expense, are probably now feasible. In pure scientific work, gas-particle DNS calculations, including the simultaneous computation of the instantaneous particle concentration and velocity fields, should help to reveal the complex secrets of particle behaviour in turbulent flows.

## Acknowledgments

The authors wish to thank Dr S. A. Slater for providing the original Lagrangian particle code. They are also grateful to Dr J. E. Fackrell of PowerGen for his technical advice. D.P.H. was supported by a Cambridge University Domestic Research Studentship with maintenance contributions from the PowerGen Power Technology Centre, Ratcliffe-on-Soar, and the Ford of Britain Trust of the Cambridge University Engineering Department.

## Footnotes

- Received April 7, 2004.
- Accepted October 1, 2004.

- © 2005 The Royal Society