Differential formulation of the viscous history force on a particle for efficient and accurate computation

It is well known that the computation of the Basset-like history force is very demanding in terms of CPU and memory requirements, since it requires the evaluation of a history integral. We use the recent rational theory of Beylkin & Monzón (Appl. Comput. Harmon. Anal., vol. 19, 2005, pp. 17–48) to approximate the history kernel in the form of exponential sums to reformulate the viscous history force in a differential form. This theory allows us to approximate the history kernel in terms of exponential sums to any desired order of accuracy. This removes the need for long-time storage of the acceleration histories of the particle and the fluid. The proposed differential form approximation is applied to compute the history force on a spherical particle in a synthetic turbulent flow and a wall-bounded turbulent channel flow. Particles of various diameters are considered, and results obtained using the present technique are in reasonable agreement with those achieved using the full history integral.


Introduction
The modelling of disperse multiphase flow often takes the form of a mixed Eulerian-Lagrangian description in which the fluid is treated as a continuum while the particles (or drops, or bubbles) are tracked by integrating an equation of motion of the form (we omit the body, buoyancy and lift forces for brevity) m * p du * p dt * = 6πa * µ * u * − u * p + m * fluid, a * the particle radius and ρ * , µ * and ν * = µ * /ρ * the density, viscosity and kinematic viscosity of the fluid. In this work, the quantities with the superscript * denotes dimensional quantities and the terms with subscript τ * indicate that the expression is to be evaluated at τ * . The terms in the right-hand side of (1.1) are the Stokes drag force (which may be corrected for finite-Reynolds-number effects), the force due to the ambient pressure gradient and the added-mass force. The last term is the so-called history force, which originates from the diffusive contribution to the transport of the vorticity generated at the particle surface. The presence of a contribution of this type was initially recognized by Boussinesq (1885) and Basset (1888), who found for the kernel K the explicit expression .
(1.2) Equation (1.1) with K = K BH is often referred to as the Basset-Boussinesq-Oseen (BBO) equation. An extension of the BBO equation which accounts for the variation of the ambient flow was given by Maxey & Riley (1983) and Gatignol (1983), who collected and systematized earlier contributions by several earlier authors starting with Basset (1888), Taylor (1928), Tchen (1947), Landau & Lifschitz (1987) and others. The BBO and Maxey-Riley-Gatignol (MRG) equations are rigorously valid only as long as the vorticity has not diffused out of the Stokes region. For later times, as recognized by Mei (1993), Lovalenti & Brady (1993a,b), the slow decay inversely proportional to √ t * − τ * turns to a faster one. For a particle starting from rest with a constant velocity, the steady Stokes drag is eventually approached as (t * − τ * ) −2 . In the case of small step increase of velocity, the transient force decays proportionally to (t * − τ * ) −5/2 e b(t * −τ * ) , with b an appropriate inverse time scale, and so on. The precise details of the late-time decay are of relatively minor importance as, by the time they set in, the kernel K has already decreased very substantially and little difference would be found in practice by using one or the other form of the kernel K, provided the early time inverse proportionality to √ t * − τ * is respected. Mei & Adrian (1992) proposed the use of a kernel of the form where f H = 0.75 + 0.105Re, Re = 2a * |u * − u * p |/ν * , and are evaluated at t * . As discussed in Lovalenti & Brady (1993a,b), at finite Reynolds number the kernel depends on both the current time t * as well as the time difference t * − τ * to account for the nonlinear inertial effect of the fluid. This dependence is reflected in the above Mei & Adrian kernel.
The kernel (1.3) reduces to the Basset form (1.2) for small times and transitions to a decay proportional to (t * − τ * ) −2 for later times. As shown by Mei & Adrian (1992), this feature captures in an adequate way the late-time increase of the decay rate of the kernel K for a large variety of situations. Whatever the exact form of K, the precise computation of the history force is expensive due to the integral representation, which destroys locality in time. The evaluation of the integral requires storage of particle and fluid acceleration histories at the particle location. Systems involving large numbers of particles will therefore incur an substantial increase in memory requirement. For this reason, the history force has often been neglected by many researchers; for example Wood & Jenkins (1973), Lee & Hsu (1994), Schmeeckle & Nelson (2003). However, there are instances where the history force has been shown to be very important; for example, Vojir & Mchaelides (1994), Nino & Garcia (1998), Mordant & Pinton (2000), Sobral, Oliveira & Cunha (2007).
There are primarily two issues involving the calculation of the history force -(i) the singularity of the Basset kernel, K BH , as (t * − τ * ) → 0 and (ii) the cost of computing the long-time history effect. Brush, Ho & Yen (1964) were the first to introduce the idea of splitting the history integral into two parts: one involving the singularity and the other free of singularity. They handled the integral involving the singularity by assuming that the acceleration was constant for a short period of time, subsequently leading to an analytical integration. A trapezoidal rule was used to evaluate the second part of the integral. Nevertheless, the cost involved in computing the second part was unresolved. The literature documents several efforts attempting to reduce the cost of the history force evaluation, focusing either on various approximations of the kernel or a simple truncation. The first approach was followed, among others, by Michaelides (1992), who reformulated the particle equation of motion (1.1) into a second-order integro-differential form explicit in the particle velocity without, however, eliminating the memory effect of the particle acceleration. Dorgan & Loth (2007) proposed to truncate the history integral after a finite time. With this approach, the accuracy of the evaluation depends on the truncation time, which, due to the slow decay of the Basset kernel, must be chosen rather large so that the cost saving is limited. Bombardelli, Gonzalez & Nino (2008) exploited the fractional derivative nature of the Basset force to obtain some small improvement in the efficiency of the computation.
A different approach was taken by van Hinsberg (2011), who proposed the use of a kernel which switches from the Basset form to a sum of exponentials. The contribution of exponential functions can be computed in a recursive manner, making evaluation of the integral efficient. Daitche (2013), using integration by parts, rewrote the history force by taking the time derivative outside the integral. Subsequently, the relative velocity was approximated using polynomials, and several higher-order quadrature schemes were presented to evaluate the resulting integral, thereby avoiding the singularity of the kernel. A somewhat improved approach to the work of Dorgan & Loth (2007) has recently been proposed by Elghannay & Tafti (2016), who expressed the history force as a product of a decaying function and an accumulated averaged acceleration. Here the decaying function depended only on the number of time steps taken. Further, the averaged acceleration is considered only up to times where the kernel behaves as 1/ √ t * − τ * (in the limit of zero Reynolds number). All of the above-mentioned formulations require the evaluation of the history integral. In the current work, however, we present a new approach which avoids explicit integration of the history kernel and thereby the requirement to store the long history of past relative accelerations. We achieve this objective by approximating the kernel, to arbitrary accuracy, in terms of exponential sums, following the recent theory of Beylkin & Monzón (2005). As will be demonstrated, there are advantages in expressing the history kernel in terms of exponential sums. One can equivalently view the approximation by exponential sums as a rational approximation of the kernel in the frequency domain. This connection will be demonstrated in § 6.
We base our approach on the Mei & Adrian (1992) kernel, which, as noted above, has the desirable feature of a late-time fast decay. The methodology presented here applies equally for any other form of the kernel. In order to accurately and efficiently capture the singularity of the kernel, we break up the integral into two parts, a 'recent past' contribution from t * − t * 0 to t * , and the remainder, which may be called the 'distant past' contribution. Efficiency is achieved by choosing a suitably small time interval t * − t * 0 involving only a few time steps. The integration over the distant past is simplified by using the exponential approximation to the kernel, which can be converted to a differential form. We establish the range of time scales of importance encountered in practical particulate turbulent flows to estimate possible values of t * 0 . We present approximations using one to four exponential terms, which results in an accuracy of up to 99 %.
By way of example, in § 2 we establish the equivalence of the history integral and the differential representation for the compressible inviscid-unsteady force. In § 3, the approximation to the viscous history force is presented using the exponential sums. The differential formulation is described in § 3.2. Section 4 discusses the range of time scales encountered in the simulation of a particle moving in a turbulent flow. We demonstrate the efficiency and accuracy of the new method with an example of synthetic turbulence in § 5. The energy implications of the approximation in the frequency space are discussed in § 6. Finally we present the conclusions in § 7.

Compressible inviscid-unsteady force
Before we look at the history force, it is instructive to consider its inviscid counterpart, the added-mass force. In an incompressible flow the added-mass force is proportional to the instantaneous relative acceleration as given in (1.1). Here, we will consider a compressible flow, where the inviscid-unsteady force assumes an integral form due to the finite speed of propagation of sound. For a quiescent ambient, this force can be written in the following non-dimensional form (Parmar, Haselbacher & Balachandar 2008, 2011Parmar, Balachandar & Haselbacher 2012a,b) in which all quantities are non-dimensionalized using a * as the length scale, a * /c * 0 as the time scale, where c * 0 is the speed of sound, and m * f as the mass scale. K IU (t) is the compressible inviscid-unsteady history kernel. For Ma → 0, where Ma = |u * p |/c * 0 is the Mach number, the compressible inviscid-unsteady history kernel can be obtained analytically (Longhorn 1952;Parmar et al. 2011) and expressed as where i = √ −1. In this limit, the compressible inviscid-unsteady kernel naturally takes the form of a sum of two exponentials with complex arguments.
Using (2.2), the Fourier domain representation of (2.1) can be written as where ω is the frequency. The transfer function connecting the particle acceleration and the compressible inviscid-unsteady force occurs as a rational function in the frequency space. Thus, the compressible inviscid-unsteady force can be exactly written as a second-order differential equation in the time domain as The integral form of compressible inviscid-unsteady force (2.1) and the differential form (2.4) are equivalent, but it is computationally advantageous to use the differential form since it does not require storage of past history of particle motion. The dependence of the force on past history is implicitly taken into account by the solution of the differential equation without the need of explicit storage. Note that to solve (2.4) we require appropriate initial conditions.
Since the kernel given in (2.2) decays exponentially, if a particle motion is unsteady on a time scale much longer than the acoustic time scale, the particle acceleration in (2.1) can be considered a constant and taken out of the integral and the history force can be written as valid for a slowly accelerating particle.

Viscous-unsteady force
In the previous section we established the equivalence of integral and differential formulations with the example of the compressible inviscid-unsteady force. In this section we will obtain a similar differential formulation for the viscous-unsteady force with the history kernel of Mei & Adrian (1992) to illustrate the method. But there is an important difference. The integral representation of the viscous history force is not a convolution integral and as a result the approach of § 2 using Fourier transform is not possible. In what follows we will first approximate the viscous history kernel as sum of decaying exponentials and follow the approach detailed in the appendix A to obtain an equivalent differential form. This approach is rigorous in the linear limit of Re → 0, and will be approximately extended to finite Re. Also, the approach to be discussed below has general applicability and can be used for other decaying history kernels.
From now on we consider non-dimensional quantities with the particle radius as length scale. We introduce the following time and mass scales where Re 0 is the Reynolds number corresponding to a reference velocity and with g H0 being evaluated at Re = Re 0 . With this non-dimensionalization, the history kernel of Mei & Adrian (1992) takes on a simplified form and the viscous-unsteady force can be written as We have written this equation for the case of an arbitrary particle motion in an otherwise quiescent ambient medium. The analysis of this section can be easily extended to the more general case where the ambient flow acceleration is non-zero.
In the above non-dimensional kernel, the factor r is a function of time since g H will vary as the particle Reynolds number Re changes over time. If the reference velocity is chosen to be the initial velocity difference between the particle and the fluid, r = 1 at the initial time t = 0. As relative velocity changes over time the factor r will vary accordingly. However, in the evaluation of the integral over ξ , r(t) is a constant. In the evaluation of (3.5) at a later time t it may be advantageous to rescale with an appropriate local reference value of relative velocity. The time and mass scales based on the local reference Reynolds number are then where g H is evaluated at t. With this local reference, the expressions for the viscousunsteady force (3.5) and the kernel (3.4) remain the same, except for the simplification r = 1.
In the limit ξ 1, the kernel (3.4) reduces to the Basset kernel ∝ ξ −1/2 and for long-time (ξ 1) the kernel decays faster as ξ −2 to account for nonlinear effects. The change from the short-time ξ −1/2 behaviour to the long-time ξ −2 behaviour occurs around ξ ∼ r −4/3 . For moderate Reynolds numbers (Re < O(10 3 )), the transition to an exponential decay mentioned earlier occurs when the kernel is already quite small. Therefore, the kernel given in (3.4) is adequate for the present discussion.
While the time integral of the Basset kernel is non-convergent, the long-time faster decay of K VU makes this kernel integrable at infinity. For a particle whose time scale of acceleration is much slower than O(T * Re ), the acceleration term in (3.5) can be taken out of the integral and K VU can be integrated to obtain the following explicit expression for the viscous-unsteady force: However, often the time scale of acceleration is of the same order as O(T * Re ) and therefore the above approximation becomes inappropriate and the integral (3.5) must be evaluated. The singularity of the kernel of (3.5) as ξ → 0 is important and must be accurately accounted for. For this purpose, as already mentioned, we split the integral into two parts, writing (3.11) The first term F VU,S contains the singularity of K VU and will be retained in the integral form. The integral can be evaluated by using low-cost techniques such as the trapezoidal rule. To keep the overall cost to a minimum, t − t 0 should be chosen equal to a few discrete-time intervals -for example, for a constant time step t, t − t 0 = n t, where n is a small integer. The acceleration history du p /dt needs to be stored only for the previous n − 1 steps. Our experience shows that n 2 is typically adequate. The focus of the present work is the treatment of the second term (3.11), to which we now turn.
3.1. Approximation of shifted kernel by exponential sums The second term in (3.9) can be written as which represents the integral of the acceleration with the shifted kernel K VU (t, ξ + t 0 ).
In what follows we will present a differential approximation for the above integral for a simpler kernel with a constant factor r = 1, and hence the explicit dependence of the kernel on t is suppressed. Generalization of this approximation to kernels with r = 1 will be considered in § 3.3. The time-shifted kernel K VU (ξ + t 0 ) for various values of t 0 = {0, 10 −2 , 10 −1 , 1, 10, 10 2 } is shown in figure 1 for r = 1. With a zero time shift, t 0 = 0, the singularity of the viscous-unsteady kernel is retained and it asymptotes to ξ −1/2 for small ξ . For t 0 > 0 the kernel K VU (ξ + t 0 ) is finite for small ξ and asymptotes to K VU (t 0 ). The removal of the singularity makes it possible to get an accurate representation of the time-shifted kernel in terms of a small number of exponentials for infinite time. Since K VU (ξ + t 0 ) as given in (3.4) asymptotes to a power law decay ξ −2 for large time, it cannot be approximated by the summation of decaying exponentials. Thus, we approximate K VU (ξ + t 0 ) on a large but finite interval, say [0, T]. We seek an approximation of the following form (3.13) In § 3.2 we will use these decaying exponentials in pairs to obtain equivalent secondorder differential forms. Therefore we have chosen an even number of terms. The procedure of Beylkin & Monzón (2005) for the evaluation of M, a k and b k can be summarized as follows. Discretize the domain [0, T] into 2N equal intervals of size T/(2N) by means of points ξ n = nT/2N with 0 n 2N. The coefficients (a k ) and the exponents (b k ) are to be chosen in such a way that, for a given M, the relation is valid at the instants ξ n for a given accuracy > 0. This objective is achieved through the following four steps: Here V denotes the complex conjugate of V. The N + 1 singular values are all real and non-negative and therefore can be sorted in a decreasing order such that σ 0 σ 1 · · · σ 2M · · · σ N . The decay rate of the singular values determines the number of exponential terms to be used in the approximation for a desired accuracy . The error is controlled by choosing M so that the singular value σ 2M , with 2M N, satisfies σ 2M < < σ 2M−1 . In order to control the error to the desired level the number of terms 2M in the exponential sum may need to be adjusted.
The roots γ k , 1 k N, of this polynomial are used to obtain the coefficients a k by solving the following least- Since the function to be approximated by the exponential sums, i.e. K VU , is realvalued, the 2M roots lie inside the unit disk. These 2M γ k values which will be used to calculate the exponents are all real and positive. The corresponding a k are all real as well. (iv) Typically, only 2M coefficients a k have absolute value larger than the desired accuracy , thus, only 2M coefficients a k need to be considered. The exponents are given by b k = −2N/T log γ k , which are all real as well.
Here we apply the above four-step procedure to the kernel K VU (ξ + t 0 ) for different values of t 0 . Singular values σ n in the approximation of K VU (ξ + t 0 ) for t 0 = 0.1 over 0 ξ T = 10 2 taking 2N + 1 = 501 equispaced samples are shown in figure 2 for various values of N and M.
For small values of t 0 the error increases, since the four-term approximation is not adequate very near the singularity. For t 0 = 10 −3 the relative error E(t 0 )/B(t 0 ) is approximately 38 %, but it rapidly falls with increasing t 0 . For t 0 = 10 −2 the relative error is only 0.5 % and falls to 0.02 % for t 0 = 0.1.
3.2. Differential formulation of F VU,L A two-term approximation, i.e. 2M = 2, will result in a differential form similar to (2.4) relating force and its first and second time derivatives to the shifted acceleration and its first derivative. In general, a 2M-term approximation will result in a differential equation relating force and its derivatives up to order 2M, to acceleration and its derivatives up to 2M − 1. Such a general expression is derived in appendix A. For 2M > 2 a higher-order differential equation is obtained. For second-order accuracy in time, second-order differential equations are optimal in terms of computational cost and storage. In particular, the form of the differential equation given in (2.4) is easier to numerically discretize. Therefore, we exploit the linearity of the problem and express the 2M-term approximation as a summation of M kernels as follows As discussed in appendix A, correspondingly the viscous-unsteady force F VU,L (t) can be written as a summation of M contributions as Note that each force contribution F k arising from the integrals of a single pair of exponential terms can be expressed as a second-order differential equation as (3.20) TABLE 2. Particle locations, particle diameters, the chosen t 0 and corresponding L 2 error norms when force is approximated using the differential form.
3.3. Generalization for time dependence In the previous two subsections we considered an approximation to the viscous kernel as sum of decaying exponentials and the resulting differential formulation of the viscous-unsteady force with the restriction r = 1. But as discussed in § 3, the value of r in the kernel (3.4) depends on the time t at which the force is evaluated. To account for this dependence we first follow the steps presented for r = 1 and obtain the exponential approximation for the kernel for different values of r. For example, table 2 shows the values of the coefficients a 1 to b 4 for different values of r at a fixed t 0 = 10 −2 . As can be seen from the table, the coefficients slowly and smoothly change with r. Based on this weak dependence we propose a simple approach to approximately account for the r-dependence of the viscous kernel. In this approach the differential form given in (3.20) will still be used for calculating the force contributions F k . However, the coefficients are taken to be time-dependent and the instantaneous value of these coefficients are obtained from the corresponding value of r. From the derivation presented in appendix A it can be verified that this approximation is tantamount to neglecting the effect of da k /dt and db k /dt in the differential form.
In summary, we propose to evaluate the history force, F VU , in two parts, as given in (3.9), that separate the contributions from the 'recent past' and from the 'distant past'. We choose to evaluate the integral (3.10) for the contribution from the 'recent past' using the acceleration history only for the last one or two time steps. Since the analytic behaviour of the kernel is known, this short-time integral can be accurately evaluated using standard integration techniques. Care must be taken to account for the singular nature of the kernel as t − ξ → 0. Evaluation of contribution from the 'distant past' amounts to solving M second-order differential equations of the form (3.20). For the case of 2M = 4, F VU,L (t) = F 1 (t) + F 2 (t), where using the values of coefficients given in The above equation is appropriate only in the limit r = 1. As r drifts away from this initial value, the value of the coefficients must accordingly be changed over time. Similar approach can be followed for other values of t 0 . A second-order accurate finite-difference approximation is adequate for solving (3.21), which requires storage of both the force and the acceleration from last three time steps. Typically an approximation with M = 2 is sufficient for accurate evaluation, as will be shown in § 5. Time integration of the above differential equations, even with the time-dependent coefficients, requires the storage of force and acceleration for only three previous time steps and can be implemented in a computationally efficient manner. Thus there can be some advantage compared to direct evaluation of the viscous-unsteady force from the history integral that requires storage of the acceleration history for many more time steps.

Time scales of interest in turbulent flows
As was shown before, the exponential approximation to K VU (ξ + t 0 ) depends on the value of t 0 , which is of the order of the time step used in the simulations and which, therefore, depends on the time scales of the flow being computed. It is therefore necessary to estimate the range of time scales encountered in particulate flows.
We consider a particle moving in a turbulent flow characterized by the Kolmogorov length and time scales l * κ and t * κ and by the integral length and time scales l * L and t * L . The ratio of the largest to the smallest time scales is t * where d * = 2a * is the particle diameter. We now use the analysis of Balachandar (2009) and Ling, Parmar & Balachandar (2013) to estimate the Reynolds number for a particle freely moving in turbulence in response to hydrodynamic forces. For simplicity of analysis, here we ignore the effect of gravity and other external forces on particle motion. If needed, their effects can be included in the analysis following Balachandar (2009). The particle time scale is where ρ * p is the density of the particle and Φ(Re) = 1 + 0.15 Re 0.687 is the finite-Reynolds-number correction to the Stokes drag (Clift, Grace & Weber 1978). Three different regimes can be identified depending on the value of particle time scale relative to t * κ and t * L (see also Balachandar & Eaton 2010), namely (I) t * p < t * κ , (II) t * κ < t * p < t * L and (III) t * L < t * p . In regime I, the particle Reynolds number is controlled by the Kolmogorov eddies. In regime II, the dominant contribution to the particle Reynolds number is controlled by eddies with a time scale of the same order as that of the particle. In regime III, the maximum Reynolds number of the particle is dictated by the velocity scale of the integral scale eddies. Following Balachandar (2009), the expressions for the particle Reynolds number in the three regimes are . Kolmogorov time scales for various density ratios (0, 2.5, 25, 1000). Lines are used for l * L /l * κ = 10 3 and square symbols for l * L /l * κ = 10 5 . Filled circle symbols are used to mark the lower limit on d * /l * κ corresponding to Stokes number of 0.1.
given by These relations, in conjunction with (4.1), show that the time-scale ratio t * κ /T * Re depends only on the particle size relative to the Kolmogorov scale, d * /l * κ , the density ratio ρ * p /ρ * and the turbulence Reynolds number Re L . Figure 4 shows the time-scale ratio as a function particle size relative to the Kolmogorov scale (i.e. versus d * /l * κ ) for various particle-to-fluid density ratios (0, 2.5, 25 and 1000) and for two values of l * L /l * κ = (10 3 , 10 5 ), or equivalently for two values of turbulence Reynolds numbers Re L = (10 4 , 10 20/3 ). We see that the dependence on Re L is very weak.
Also of relevance is the particle Stokes number St κ = t * p /t * κ , defined in terms of the Kolmogorov scale. Since the Kolmogorov time scale is smaller than the time scale of all other turbulent eddies, St κ is guaranteed to be larger than any other definition of Stokes number based on other turbulent scales. In other words, if St κ 1 then particle time scale is much smaller than all turbulence time scales. As discussed in Balachandar (2009) and Ling et al. (2013), St κ is also a function of d * /l * κ and particleto-fluid density ratio and again weakly dependent on Re L .
Here we are interested in evaluating the viscous history force only on particles whose the Stokes number is larger than a small threshold, say St κ > 0.1. For computing the motion of smaller particles whose Stokes number is small one need not solve the integro-differential (1.1) or the MRG equation. In the limit of St κ → 0 the particle instantly responds to all turbulent scales and the particle behaves as a tracer with the particle velocity the same as the local fluid velocity. For non-zero Stokes number the particle velocity deviates from the local fluid velocity. However, for small Stokes number (e.g. St κ < 0.1) the particle velocity can still be explicitly expressed in terms of the local fluid velocity and its space and time derivatives (Ferry & Balachandar 2001, 2002Ferry, Rani & Balachandar 2003). This equilibrium Eulerian solution for the particle velocity is in fact the asymptotic solution of the BBO equation. Thus, explicit evaluation of the viscous history force and its inclusion in (1.1) for tracking the motion of a particle becomes important only for particles whose St κ > 0.1.
In figure 4 we mark on each curve the point where St κ = 0.1 with a filled black circle. These correspond to the lower limits of d * /l * κ ; and are obtained by setting Only the portion of time-scale ratio versus (d * /l * κ ) for which St κ > 0.1 is plotted. For particles smaller than the cutoff value of d * /l * κ , indicated by the filled black circle, the Stokes number is so small that viscous history force is not of importance. Thus figure 4 focuses on the parametric space where evaluation of the viscous history force is of importance. From the figure it is clear the time-scale ratio t * κ /T * Re ranges from 10 −2 to 10 2 . In direct numerical simulations of turbulent flows the time step is chosen to be of the order of Kolmogorov time scale. Thus, over the entire parametric range within which evaluation of the viscous history force can be of importance, the scaled time step ( t * /T * Re ) ranges only from 10 −2 to 10 2 . Since we restrict t 0 to be only one or two non-dimensional time steps, in the regime of interest, the value of t 0 can be expected to be also in the range 10 −2 to 10 2 . As shown in table 1 and figure 1, for this range of t 0 value the viscous history force is very well approximated by the differential form.

Particle subjected to synthetic turbulence
We demonstrate the accuracy of approximating of the exact history integral with the present simplified differential approximation in a synthetic turbulent flow. The velocity field seen by the particle is approximated as a superposition of non-dimensional frequencies ω j from j = 1 to 100 which can be written as where ζ j and θ j are randomly generated coefficients. Figure 5 compares the viscousunsteady force obtained using the integral (3.5) and the 2M = 4 differential formulation ((3.21) corresponding to t 0 = 10 −2 ) demonstrating good agreement. Note that (3.10) was used to evaluate the short-time contribution, which remained the same for both the integral and differential formulation. The simplified formulation slightly overpredicts the viscous-unsteady force with the L 2 -norm of the error at 3 %. But the simplified approach is memory efficient and computationally faster. . (Colour online) Computation of the history force using the integral form (3.5) and the differential form (3.21) using a four-term approximation corresponding to t 0 = 10 −2 in table 1 for a particle moving with velocity given by (5.1).

Particle fixed in a turbulent channel flow
Here we consider the history (viscous-unsteady) force experienced by a spherical particle situated in a turbulent channel flow. The Reynolds number based on the half-channel height and friction velocity is Re τ = 400. The velocity fluctuations seen by the particle are taken from the work of Lee, Ha & Balachandar (2012). The non-dimensional particle location (y + ) and particle size (d + ) considered in this study are summarized in table 2; and in all cases the particle is held fixed in position. Thus we consider the case of a wall turbulence sweeping over a particle. Here, y + = y * u * /ν * and d + = d * u * /ν * , where u * is the friction velocity based on the time-averaged wall shear stress. We note that the presence of the channel walls requires d + 2y + , based on which the values in table 2 have been chosen. For the four cases chosen in this study, the history force is evaluated using the integral and differential forms presented in § § 2 and 3, respectively. However, the force is rescaled and given by where Re = d * u * p /ν * is the time-averaged particle Reynolds number. We begin by evaluating the force on a particle with a diameter, d + = 7.8 situated at y + = 3.9. The history force computed using (3.5) is shown in figure 6(a); the force is fluctuating owing to the turbulent flow velocity. Also shown in the figure is the approximate history force computed with the integral where r is held fixed at its initial value of r = 1. In this case, as the turbulent flow varies over the particle the value of r ranges from 0.5 to 1.5 and the corresponding effect of time-varying FIGURE 6. (Colour online) (a) Total viscous-unsteady force, (b) the early and late-time contributions to the history force and (c) comparison of the integral and differential forms (two-term and four-term) used in evaluating the late-time force contribution for the case of y + = 3.9 and d + = 7.8.
r is small. The r.m.s difference between the results for time-varying r and r = 1 approximation is only 1.5 %. Similarly the difference between the constant coefficient and time-varying coefficient of the differential formulation is not significant. In order that we clearly observe the contribution to the total force from recent past and distant past, we focus on some finite time range, for example, 5500 t 6000 and the results shown in figure 6(b). For evaluating the contribution from relative acceleration of the recent past, we use (3.10), where the singularity occurs at the current time.
Here we have chosen t 0 = t, which happens to be 0.17 (DNS data, Lee et al. 2012). Further, the distant past contribution is computed using (3.21); however, with two modifications: (i) coefficients of the differential equation are replaced with those corresponding to t 0 = 0.17 and (ii) particle velocity is replaced by fluid velocity. The resulting behaviour is shown in figure 6(b), where one observes both the recent past and distant past contributions are equally important. In figure 6(c) we show how the distant past contribution to the history integral compares with the differential equation approximation. If a two-term (2M = 2) approximation is used we observe that the prediction is poor, with an L 2 error norm of 41 %. However, the four-term (2M = 4) approximation leads to a much improved prediction with L 2 error norm dropping to just 8 %. It must be mentioned that the data in table 1 is used for interpolating the coefficients a 1 through b 4 for any value of t 0 that lies between 10 −3 and 10. Note that here L 2 error norm is computed as where the subscripts Int and Diff denote that the forces are evaluated using the integral and differential approximations, respectively, · 2 represents the L 2 -norm. The history force experienced by a particle situated at y + = 10.25 with d + = 1 is shown in figure 7(a). As can be seen from the figure, the force contribution from the recent past (root mean square of F + VU,S = 0.631) dominates the distant past contribution TABLE 3. Time taken (per time step, per particle, per space dimension) in milliseconds, to evaluate the history force on a particle situated in a turbulent channel flow on using the full history integral and the differential form equivalent (M = 1 and M = 2). .
(root mean square of F + VU,L = 0.0631). As was shown in figure 1, the history effects can be considered negligible for times t 3. For the case under consideration (y + = 10.25, d + = 1), we see from table 2 that t 0 = t = 5, and as a consequence the contribution from the distant past (contribution to the history integral for times −∞ t t − t) is rendered irrelevant. Since the time scale of particle acceleration is much larger than T * Re , one may choose to evaluate the force using (3.8) and the results shown in figure 7(b). As can be seen from the figure, the force evaluated using the total history integral is nearly identical to that achieved using the limiting value (3.8). Now if we consider a particle at the same location -however, with a larger diameter (y + = 10.25, d + = 10.25) -both force components (F + VU,S , F + VU,L ) become equally important since t 0 = t reduces to 0.14. The distant past contributions computed using both the integral and differential formulations are plotted in figure 7(c). As can be seen from the figure, the prediction using the four-term differential form approximation captures the force well with an L 2 error norm of 6 %. Similar behaviour is also observed for the case of y + = 400 with d + = 10, as shown in figure 7(d). In order to quantify the computational cost of using the differential form presented in this article, we compute the time taken to evaluate the long-time history force (F + VU,L ) by solving the full history integral and compare it against the differential form equivalent (3.20). We present the results in table 3 for both M = 1 and M = 2 (corresponding to the one-term and two-term approximations, respectively) and compare against the integral form. The times shown are evaluated per time step, per particle, and per space dimension. As can be seen from table 3, using the current formulation reduces the time required for evaluation for all choices of y + and d + considered. Note that the case y + = 10.25, d + = 1 is not shown, as the long-time history force was found to be negligible compared to the short-time history force (F + VU,S ). It must be cautioned that computational cost of the differential form will increase with the inclusion of more terms in the expansion, while the cost of the integral form can be decreased by limiting the integration range.

Energy implication
If we ignore the nonlinear effect and assume the history force to be a convolution integral, the exponential sum approximation has a transparent interpretation in the frequency domain where (3.12) becomes where F [(·)] denotes the Fourier transform. The transform of the history kernel appears as a multiplicative transfer function that relates the acceleration to the history force. With the approximation (3.13) we have which shows that the exponential approximation is equivalent to a rational approximation of F [K VU (t + t 0 )] in the frequency space with numerator and denominator polynomials of degrees 2M − 1 and 2M, respectively. Expressing (6.1) as we readily obtain the differential equation shown in appendix A. In fact, one could attempt to approximate directly the kernel in the frequency domain by a rational function and evaluate the optimal coefficients of the polynomials in the numerator and the denominator. This procedure, however, would involve approximating the kernel in the complex plane, which is a more difficult proposition than on the real line. For this reason, the formal theory of Beylkin & Monzón (2005) is preferable.
In the context of a point-particle approach, the change in kinetic energy of a particle due to the quasi-steady drag force F qs is given by F qs · u p . The corresponding contribution to the change of the fluid kinetic energy is −F qs · u. Although the force on the particle and the surrounding fluid have an opposite sign, the corresponding contributions to the kinetic energy of the particle and fluid do not add to zero. The net contribution to rate of kinetic energy, F qs · (u p − u) is negative definite. This dissipation of kinetic energy contributes to the fluid internal energy. An analogous consideration, in the case of the viscous history, force is not straightforward in the time domain. But the analysis is greatly facilitated by shifting the argument to the Fourier domain.
It is instructive to consider first the compressible inviscid-unsteady kernel (2.2), whose Fourier transform is Due to causality the real and imaginary parts are related by the Kramers-Kronig relation as (see page 112 of Prosperetti (2011)): where Re{} and Im{} denote real and imaginary parts of a complex quantity and PV denotes the principal value. When multiplied by the particle velocity, the real part of the inviscid-unsteady force gives rise to a reversible change of the kinetic energy of the fluid. In the limit of incompressible flow, ω → 0 and correspondingly F[K IU ] → 1/2. Thus, in this limit, F[K IU ] is purely real and the work done on the fluid is entirely reversible. In a compressible flow, with the introduction of a finite speed of sound, the imaginary part of F[K IU ] becomes non-zero, and the corresponding contribution to the force accounts for the irreversible loss of energy by acoustic radiation. In a similar manner, the Fourier transform of each K k approximation of the viscousunsteady kernel (see (3.18)) is Here again the real part corresponds to the reversible component and arises from the viscous modification of the near-field, while the imaginary part relates to the added viscous dissipation due to the viscous-unsteady force. The separation of reversible and irreversible components of the history force is not obvious in the time domain. For example, if (2.1) and (3.5) are used to compute inviscid-and viscous-unsteady forces as a function of time, at any given instance the computed force cannot be separated into reversible and dissipative contributions. This difficulty arises from the fact that the frequency decomposition of relative acceleration at any given instant depends on how it changes in the future. Thus, reversible and dissipative contributions of F IU (t) and F VU (t) can be determined after times of O(a * /c * 0 ) and O(T * Re ) when the respective kernels decay sufficiently such that they do not contribute to the history integral. As discussed in § 2, for a particle accelerating on a time scale much slower than O(a * /c * 0 ), the compressible inviscid-unsteady force reduces to the added-mass force. All the work done by the compressible inviscid-unsteady force goes to change the kinetic energy of the near-field. Similarly, for a particle accelerating on a much slower time scale compared to O(T * Re ), the transform of the history kernel reduces to 8π/(9 √ 3) and the work done by the viscous-unsteady force goes towards the reversible kinetic energy content of the near-field. Viscous dissipation will be entirely accounted for by the slowly varying quasi-steady force.

Conclusion
By expressing the history kernel as exponential sums we have developed an approximate differential formulation of the viscous-unsteady force which provides an alternate approach to the traditional history integral. We first illustrate this with the compressible inviscid-unsteady force, where the kernel naturally occurs as a sum of two exponentials. As a result the history integral and the differential formulation of the compressible inviscid-unsteady force are identical. Recent rational theory of Beylkin & Monzón (2005) allows us to approximate the viscous history kernel in terms of exponential sums to any desired order of accuracy. To overcome difficulties due to the singularity of the viscous history kernel at t = 0, the history integral is split into two parts. The first part is integration over a few time steps to capture the short-time singular behaviour of the viscous kernel accurately. The second part is expressed as a differential form using exponential sums to approximate the kernel. A four-term approximation has been shown to approximate the finite-Reynolds-number viscous-unsteady kernel of Mei & Adrian (1992) (3.4) to within 1 % accuracy. From a simple analysis of particle motion in turbulent flows, we argue that the proposed approach is appropriate for a wide range of particle properties and turbulence Reynolds numbers. The proposed differential form approximation was used to evaluate the history force on a particle in a synthetic turbulent flow and a turbulent channel flow. Particle of varying diameters located at different locations within the channel were considered. In all the cases, we found the four-term differential form to capture the viscous history adequately (L 2 error norm 6 %) and the error can be further reduced with inclusion of additional terms in the expansion. The exponential sum approximation of the history kernel makes it possible to characterize the reversible and irreversible contributions to the near-field energy.