2.1 The global relation and some notation

Next we integrate equation (2.2) over the entirety of $\mathscr{D}$ and employ the divergence theorem to obtain a ‘global relation’,

(2.3) $$\begin{eqnarray}\displaystyle & & \displaystyle \iint _{\mathscr{D}}[(\text{e}^{-\text{i}kx+\unicode[STIX]{x1D714}(k)t}\boldsymbol{q})_{t}-(\text{e}^{-\text{i}kx+\unicode[STIX]{x1D714}(k)t}(\boldsymbol{q}_{x}+\text{i}k\boldsymbol{q}))_{x}]\,\text{d}x\,\text{d}t=0,\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle \Rightarrow & & \displaystyle \int _{\text{d}\mathscr{D}}[\text{e}^{-\text{i}kx+\unicode[STIX]{x1D714}(k)t}\boldsymbol{q}\,\text{d}x+\text{e}^{-\text{i}kx+\unicode[STIX]{x1D714}(k)t}(\boldsymbol{q}_{x}+\text{i}k\boldsymbol{q})\,\text{d}t]=0,\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle \Rightarrow & & \displaystyle \int _{0}^{\infty }\text{e}^{-\text{i}kx}\boldsymbol{q}_{0}(x)\,\text{d}x-\int _{0}^{\infty }\text{e}^{-\text{i}kx+\unicode[STIX]{x1D714}(k)T}\boldsymbol{q}(x,T)\,\text{d}x-\int _{0}^{T}\text{e}^{\unicode[STIX]{x1D714}(k)t}(\boldsymbol{q}_{x}(0,t)+\text{i}k\boldsymbol{g}_{0}(t))\,\text{d}t=0.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The terms at the boundary of $\mathscr{D}$ at infinity vanish since $\boldsymbol{q}(x,t)\rightarrow 0$ as $x\rightarrow \infty$ . Let us define the Fourier transform in pseudo-space as

(2.6) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{q}}_{0}(k) & = & \displaystyle \int _{0}^{\infty }\text{e}^{-\text{i}kx}\boldsymbol{q}_{0}(x)\,\text{d}x,\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{q}}(k,T) & = & \displaystyle \int _{0}^{\infty }\text{e}^{-\text{i}kx}\boldsymbol{q}(x,T)\,\text{d}x,\end{eqnarray}$$

and the time transform of the boundary terms as

(2.8) $$\begin{eqnarray}\displaystyle \tilde{\boldsymbol{g}}_{0}(\unicode[STIX]{x1D714},T) & = & \displaystyle \int _{0}^{T}\text{e}^{\unicode[STIX]{x1D714}(k)t}\boldsymbol{g}_{0}(t)\,\text{d}t,\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle \tilde{\boldsymbol{g}}_{1}(\unicode[STIX]{x1D714},T) & = & \displaystyle \int _{0}^{T}\text{e}^{\unicode[STIX]{x1D714}(k)t}\boldsymbol{q}_{x}(0,t)\,\text{d}t.\end{eqnarray}$$

Further we define ${\mathcal{D}}^{+}$ as

(2.10) $$\begin{eqnarray}\displaystyle {\mathcal{D}}^{+}=\{k\in \mathbb{C}:\text{Re}(\unicode[STIX]{x1D714}(k))<0,~\text{Im}(k)>0\}, & & \displaystyle\end{eqnarray}$$

and ${\mathcal{D}}^{-}$ as

(2.11) $$\begin{eqnarray}\displaystyle {\mathcal{D}}^{+}=\{k\in \mathbb{C}:\text{Re}(\unicode[STIX]{x1D714}(k))<0,~\text{Im}(k)<0\}, & & \displaystyle\end{eqnarray}$$

as shown in figure 2 $(b)$ . Note that $\text{e}^{\unicode[STIX]{x1D714}(k)t}$ is analytic, bounded and decaying for large $k$ in these domains, a property we shall use repeatedly. Finally we arrive at the global relation given by

(2.12) $$\begin{eqnarray}\displaystyle \hat{\boldsymbol{q}}_{0}(k)-\text{e}^{\unicode[STIX]{x1D714}(k)T}\hat{\boldsymbol{q}}(k,T)-\tilde{\boldsymbol{g}}_{1}-\text{i}k\tilde{\boldsymbol{g}}_{0}=0,\quad k\in \mathbb{C}^{-}. & & \displaystyle\end{eqnarray}$$

Since the initial condition is zero everywhere on the real line, we have

(2.13) $$\begin{eqnarray}\displaystyle \text{e}^{\unicode[STIX]{x1D714}(k)T}\hat{\boldsymbol{q}}(k,T)+\tilde{\boldsymbol{g}}_{1}+\text{i}k\tilde{\boldsymbol{g}}_{0}=0,\quad k\in \mathbb{C}^{-}. & & \displaystyle\end{eqnarray}$$

2.2 Dirichlet-to-Neumann map

Equation (2.1) describes the Dirichlet boundary-value problem for the diffusion equation. We now compute the associated Neumann condition, i.e. we ask what Neumann boundary condition would give the same solution as this Dirichlet condition in the entire domain. To do so, we multiply equation (2.13) by $\text{i}k\text{e}^{-\unicode[STIX]{x1D714}(k)t}$ for $0<t<T$ and integrate over the contour $\unicode[STIX]{x2202}{\mathcal{D}}^{-}$ , the boundary of ${\mathcal{D}}^{-}$ shown in figure 2(b):

(2.14) $$\begin{eqnarray}\displaystyle \int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\left[\text{i}k\int _{0}^{T}\text{e}^{\unicode[STIX]{x1D714}(s-t)}\boldsymbol{q}_{x}(0,s)\,\text{d}s-k^{2}\int _{0}^{T}\text{e}^{\unicode[STIX]{x1D714}(s-t)}\boldsymbol{q}(0,s)\,\text{d}s+\text{i}k\text{e}^{\unicode[STIX]{x1D714}(k)(T-t)}\boldsymbol{q}(k,T)\right]\text{d}k=0. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

After the manipulations detailed in appendix A, we obtain

(2.15) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x03C0}\boldsymbol{q}_{x}(0,t)-\sqrt{\frac{\unicode[STIX]{x03C0}}{t}}\,\boldsymbol{q}(0,0)-\int _{-\infty }^{\infty }\int _{0}^{t}\text{e}^{-k^{2}(t-s)}\dot{\boldsymbol{q}}(0,s)\,\text{d}s\,\text{d}k=0. & & \displaystyle\end{eqnarray}$$

Setting $k^{2}(t-s)=m^{2}$ , where $m$ is real, and evaluating the resulting Gaussian integral, we obtain the Dirichlet–Neumann map as

(2.16) $$\begin{eqnarray}\displaystyle \boldsymbol{q}_{x}(0,t)=-\sqrt{\frac{1}{\unicode[STIX]{x03C0}t}}\,\boldsymbol{q}(0,0)-\frac{1}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{t}\frac{\dot{\boldsymbol{q}}(0,s)}{\sqrt{t-s}}\,\text{d}s. & & \displaystyle\end{eqnarray}$$

The integral term above is the definition of the Riemann–Liouville half-derivative. This expression relating the Neumann condition for the diffusion equation and the Riemann–Liouville half-derivative forms the basis of our reformulation of the MR equations. More general connections of boundary-value problems to partial differential equations and fractional derivatives are being explored in the authors’ work currently underway.

3.1 Solution for time-dependent forcing

The form of the MR equation (3.2)–(3.3) suggests that the forcing function is both space- and time-dependent. However, it is instructive and useful to first study a case when the forcing is a function of time alone, i.e.  $\boldsymbol{f}(t)$ . The ideas presented here will be extended to more complicated scenarios in later sections.

Using the definitions (2.8)–(2.9), we take the time transform of the boundary condition (3.6) and after integrating by parts we get

(3.10) $$\begin{eqnarray}\displaystyle (\text{e}^{\unicode[STIX]{x1D714}(k)T}\boldsymbol{q}(0,T)-\boldsymbol{v}_{0})+(\unicode[STIX]{x1D6FC}-k^{2})\tilde{\boldsymbol{g}}_{0}-\unicode[STIX]{x1D6FE}\tilde{\boldsymbol{g}}_{1}=\,\tilde{\boldsymbol{f}}. & & \displaystyle\end{eqnarray}$$

Here $\tilde{\boldsymbol{f}}$ is the time transform of the function $\boldsymbol{f}$ . Using the global relation (2.13), we eliminate $\tilde{\boldsymbol{g}}_{1}$ from the above expression and multiply the resulting equality by $\text{i}k\text{e}^{-\unicode[STIX]{x1D714}(k)t}$ for $0<t<T$ and finally integrate over $\unicode[STIX]{x2202}{\mathcal{D}}^{-}$ . This leads to

(3.11) $$\begin{eqnarray}\displaystyle \int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\text{i}k\tilde{\boldsymbol{g}}_{0}(\unicode[STIX]{x1D714},T)\text{e}^{-\unicode[STIX]{x1D714}(k)t}\,\text{d}k=\int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\text{i}k\text{e}^{-\unicode[STIX]{x1D714}(k)t}\left[\frac{\,\tilde{\boldsymbol{f}}-\text{e}^{\unicode[STIX]{x1D714}T}\boldsymbol{q}(0,T)+\boldsymbol{v}_{0}-\unicode[STIX]{x1D6FE}\text{e}^{\unicode[STIX]{x1D714}T}\hat{\boldsymbol{q}}(k,T)}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\right]\text{d}k, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\text{i}k\text{e}^{-k^{2}t}\int _{0}^{T}\text{e}^{k^{2}s}\boldsymbol{q}(0,s)\,\text{d}s\,\text{d}k\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\frac{\text{i}k\text{e}^{-k^{2}t}(\,\tilde{\boldsymbol{f}}(k,T)+\boldsymbol{v}_{0})}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}k-\int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\text{i}k\frac{\text{e}^{k^{2}(T-t)}(\boldsymbol{q}(0,T)+\unicode[STIX]{x1D6FE}\hat{\boldsymbol{q}}(k,T))}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}k,\end{eqnarray}$$

where (in the second line) we have used the definition of $\tilde{\boldsymbol{g}}_{0}$ and the fact that $\unicode[STIX]{x1D714}(k)=k^{2}$ . We substitute $k^{2}=\text{i}l$ in the left-hand-side term above, which describes the boundary $\unicode[STIX]{x2202}{\mathcal{D}}^{-}$ , see figure 2(b), to arrive at

(3.13) $$\begin{eqnarray}\displaystyle & & \displaystyle -\frac{1}{2}\int _{-\infty }^{\infty }\text{e}^{-\text{i}l(t-s)}\int _{0}^{T}\boldsymbol{q}(0,s)\,\text{d}l\,\text{d}s\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\frac{\text{i}k\text{e}^{-k^{2}t}(\,\tilde{\boldsymbol{f}}(k,T)+\boldsymbol{v}_{0})}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}k-\int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\text{i}k\frac{\text{e}^{k^{2}(T-t)}(\boldsymbol{q}(0,T)+\unicode[STIX]{x1D6FE}\hat{\boldsymbol{q}}(k,T))}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}k,\end{eqnarray}$$

(3.14) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x03C0}\boldsymbol{q}(0,t) & = & \displaystyle \int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\int _{0}^{t}\frac{\text{i}k\text{e}^{-k^{2}(t-s)}\boldsymbol{f}(s)}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}k\,\text{d}s+\int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\frac{\boldsymbol{v}_{0}\text{i}k\text{e}^{-k^{2}t}}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}k\nonumber\\ \displaystyle & & \displaystyle +\int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\int _{t}^{T}\frac{\text{i}k\text{e}^{k^{2}(s-t)}\boldsymbol{f}(s)}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}k\,\text{d}s\nonumber\\ \displaystyle & & \displaystyle -\int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\text{i}k\frac{\text{e}^{k^{2}(T-t)}(\boldsymbol{q}(0,T)+\unicode[STIX]{x1D6FE}\hat{\boldsymbol{q}}(k,T))}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}k,\end{eqnarray}$$

where we have used the standard Fourier inversion formula on the left-hand side. There are two poles due to the quadratic in the denominator of the above equation which are located at

(3.15) $$\begin{eqnarray}\displaystyle k_{\pm }={\textstyle \frac{1}{2}}(\text{i}\unicode[STIX]{x1D6FE}\pm \sqrt{4\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FE}^{2}}\,). & & \displaystyle\end{eqnarray}$$

Since both $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FC}$ are positive parameters, the poles $k_{\pm }$ always lie in $\mathbb{C}^{+}$ and reach the real line for $\unicode[STIX]{x1D6FE}=0$ , i.e. when Basset history is absent (see figure 3). An appeal to Jordan’s lemma (Ablowitz & Fokas Reference Ablowitz and Fokas2003, p. 222) shows that the third and fourth terms (those involving function evaluations at $T$ ) in equation (3.14) do not contribute, since $k^{2}$ is negative in ${\mathcal{D}}^{-}$ and $(t-s)>0$ . It is not surprising that these two terms do not contribute, as this only ensures causality: the dynamics of the particle for all times $t<T$ cannot depend on a quantity from future time $T$ . Using the fact that $\text{e}^{-k^{2}t}$ is decaying in the region outside ${\mathcal{D}}^{-}$ , we can write the final solution expression as

(3.16) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x03C0}\boldsymbol{q}(0,t)=\int _{\unicode[STIX]{x2202}{\mathcal{D}}^{-}}\int _{0}^{t}\frac{\text{i}k\text{e}^{-k^{2}(t-s)}\boldsymbol{f}(s)}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}s\,\text{d}k-\int _{-\infty }^{\infty }\frac{\boldsymbol{v}_{0}\text{i}k\text{e}^{-k^{2}t}}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}k. & & \displaystyle\end{eqnarray}$$

Note that the particle velocity is now given by the above expression directly (and explicitly) in terms of the initial condition $\boldsymbol{v}_{0}~(=\boldsymbol{v}(0)-\boldsymbol{u}(0))$ and the forcing $\boldsymbol{f}(t)$ in the relative frame of reference. The quantity $\boldsymbol{q}(0,t)$ represents the Dirichlet boundary condition for the modified Robin boundary-value problem for the diffusion equation. Once the relative velocity $\boldsymbol{q}(0,t)$ is available, one may readily compute particle trajectories $\boldsymbol{y}(t)$ from (3.7).

Figure 3.

The region in $k\in \mathbb{C}^{+}$ where the poles $k_{\pm }$ from (3.15) are located for Stokes number $S$ in the range 0.01–1 and density ratios corresponding to $R$ in the range $1/3$ –5. As the Stokes number increases, the poles approach the real line shown using arrows, and similarly as we increase the density ratio. Note that, without Basset history, the poles would merely lie symmetrically on the $k_{r}$ line.

Without the unified transform method, one may not have suspected such a boundary-value problem is indeed solvable. Of course, one could in principle have employed the Laplace transform to obtain the solution, though a couple of issues arise then. If one were to use the Laplace transform, the inverse would involve branch cuts in the general case, which are difficult to evaluate, and, as mentioned above, explicit solutions are most often not possible. However, the unified transform method provides the appropriate reparametrisation to avoid these complications, and to further provide an explicit expression in $t$ for the solution, given by (3.16). As we shall see in the following sections, especially when considering the nonlinear problem $\boldsymbol{f}(\boldsymbol{q}(0,t),\boldsymbol{y},t)$ , the true benefit of adopting the new formulation is the local nature of the extended system, which is precisely due to the connection between the MR equations and the diffusion equation.

4.1 Example 1: relaxing particle

In the absence of a background flow (namely when $\boldsymbol{f}=0$ ), expression (3.16) for the relative velocity $\boldsymbol{q}(0,t)$ may be analysed using Watson’s lemma. Since the real line in $k$ happens to be the path of steepest descent, Laplace’s method provides the asymptotic behaviour for functions of the form

(4.1) $$\begin{eqnarray}\displaystyle \mathscr{K}(\unicode[STIX]{x1D701}):=\int _{a}^{b}\text{e}^{g(\unicode[STIX]{x1D701})t}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701}\quad \text{as}~t\rightarrow \infty ,~t>0, & & \displaystyle\end{eqnarray}$$

in terms of the local maximum value of the function $g(\unicode[STIX]{x1D701})$ when $\unicode[STIX]{x1D701}\in [a,b]$ . This leads to the well-known asymptotic result for integrals of the above form:

(4.2) $$\begin{eqnarray}\displaystyle \mathscr{K}(\unicode[STIX]{x1D701}):=\int _{-\infty }^{\infty }\text{e}^{-\unicode[STIX]{x1D701}^{2}t}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D701})\,\text{d}\unicode[STIX]{x1D701}~\approx ~\sqrt{\frac{\unicode[STIX]{x03C0}}{t}}\mathop{\sum }_{n=0}^{\infty }\frac{\unicode[STIX]{x1D713}_{k}(0)}{2^{2n}n!}\,t^{-n},\quad \text{as}~t\rightarrow \infty . & & \displaystyle\end{eqnarray}$$

On setting $\boldsymbol{f}=0$ we see that the dominant contribution of the integral in (3.16) comes from $k=0$ . To leading order we have

(4.3) $$\begin{eqnarray}\displaystyle \boldsymbol{q}(0,t)\approx \frac{\boldsymbol{v}_{0}}{\sqrt{\unicode[STIX]{x03C0}}t^{3/2}}\,\frac{\unicode[STIX]{x1D6FE}}{2\unicode[STIX]{x1D6FC}^{2}}+O(t^{-5/2}),\quad t\gg 0. & & \displaystyle\end{eqnarray}$$

Farazmand & Haller (Reference Farazmand and Haller2015) and Langlois et al. (Reference Langlois, Farazmand and Haller2015) approached this problem by using a Laplace transform. Their complete solution was provided in terms of convolutions and not directly as we obtain in (3.16). For large times, however, Langlois et al. (Reference Langlois, Farazmand and Haller2015) present the same equation as (4.3). Their solution for a relaxing particle is the only earlier closed-form solution that we know of, which accounts for the initial conditions correctly. We contrast the complete solution (3.16) with the one obtained upon neglecting the Basset history integral in the original MR equation. This is equivalent to neglecting the Neumann term $\boldsymbol{q}_{x}(0,t)$ in (3.10) and repeating the calculation of $\boldsymbol{q}(0,t)$ . Alternatively, we just set $\unicode[STIX]{x1D6FE}=0$ in (3.16) to obtain

(4.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x03C0}\boldsymbol{q}(0,t)=\int _{-\infty }^{\infty }\frac{\boldsymbol{v}_{0}\text{i}k\text{e}^{-k^{2}t}}{k^{2}-\unicode[STIX]{x1D6FC}}\,\text{d}k. & & \displaystyle\end{eqnarray}$$

The non-zero contribution to the integrand comes only from the poles as the integrand is odd, leading to

(4.5) $$\begin{eqnarray}\displaystyle \boldsymbol{q}(0,t)=\boldsymbol{v}_{0}\text{e}^{-\unicode[STIX]{x1D6FC}t}, & & \displaystyle\end{eqnarray}$$

and is precisely the expression one gets by solving the MR equation directly without the history term (and $\boldsymbol{f}=0$ ).

Figure 4.

Particle velocity with (solid lines) and without (dashed lines) Basset history, along with the asymptotic solution (dotted lines) from (4.3). The curves from bottom to top, respectively, correspond to three different density ratios, at a fixed Stokes number of $S=1$ : light particles, with $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D71A}_{p}/\unicode[STIX]{x1D71A}_{f}=0.01$ ; neutrally buoyant particles, with $\unicode[STIX]{x1D6FD}=1$ ; and heavy particles, with $\unicode[STIX]{x1D6FD}=5$ . The decay rate of particle velocity goes as $t^{-3/2}$ for long times when Basset history is accounted for, in contrast to the exponential decay when it is neglected.

In figure 4 we plot the velocity obtained both by including the history integral and by excluding it for three cases: neutrally buoyant ( $\unicode[STIX]{x1D6FD}=1$ ), light ( $\unicode[STIX]{x1D6FD}=0.01$ ) and heavy ( $\unicode[STIX]{x1D6FD}=5$ ) particles. Two important features are evident. First, the solution with Basset history settles not exponentially but as a power law (in this case a very slow $t^{-3/2}$ ). Secondly, there exists a cross-over time prior to which a particle evolving under the influence of the Basset history integral and Stokes drag relaxes faster than a particle evolving only under Stokes drag. However, after the cross-over time the particle that evolves only under Stokes drag relaxes faster. This contradicts the popular notion that Basset history acts as an effective drag. Thus far, for this example, we have considered the effect of varying density ratios on the particle relaxation time. Note that varying the Stokes number does not lead to substantially different velocity profiles since the Stokes number may be scaled out of the solution by replacing $k\rightarrow \tilde{k}/\sqrt{S}$ and $t\rightarrow \tilde{t}S$ . Particle velocities for Stokes number $S\neq 1$ may be obtained from those depicted in figure 4 by suitably scaling time.

4.2 Example 2: sedimenting particle

The dynamics of particles settling under gravity is relevant in industrial applications, for aerosol particles and droplets in the Earth’s atmosphere, and for carbon sequestration by marine snow in the oceans. In this section a two-dimensional (2-D) flow is considered, with a superscript $(\cdot )^{(2)}$ indicating the vertical component of the relevant vector, while $(\cdot )^{(1)}$ indicates the horizontal component. When the particle is acted upon by a body force such as gravity, the vertical component of the forcing function $\boldsymbol{f}$ , $f^{(2)}(t)$ , is a constant of magnitude $\unicode[STIX]{x1D70E}$ . This allows us to evaluate the time transform explicitly as

(4.6) $$\begin{eqnarray}\displaystyle \tilde{f}^{(2)}(\unicode[STIX]{x1D714},T)=\int _{0}^{T}\unicode[STIX]{x1D70E}\text{e}^{\unicode[STIX]{x1D714}s}\,\text{d}s=\frac{\unicode[STIX]{x1D70E}(\text{e}^{\unicode[STIX]{x1D714}T}-1)}{\unicode[STIX]{x1D714}}. & & \displaystyle\end{eqnarray}$$

The settling velocity of the particle is

(4.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x03C0}q^{(2)}(0,t)=\int _{-\infty }^{\infty }\frac{\text{i}k\text{e}^{-k^{2}t}v_{0}^{(2)}}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}k+\int _{-\infty }^{\infty }\frac{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FE}(1-\text{e}^{-k^{2}t})}{(\unicode[STIX]{x1D6FC}-k^{2})^{2}+k^{2}\unicode[STIX]{x1D6FE}^{2}}\,\text{d}k. & & \displaystyle\end{eqnarray}$$

A solution for this problem is provided in Clift et al. (Reference Clift, Grace and Weber2005). However, their solution necessitates a condition of zero initial velocity for the particle. Moreover, their solution is in terms of sums of error functions, and so the asymptotic form for large time is difficult to glean directly. We evaluate the full solution (4.7), without the requirement for zero initial velocity, for the same particle density ratios as in the previous example (but fixed Stokes number $S=0.01$ ), corresponding to neutrally buoyant, light and heavy particles: $\unicode[STIX]{x1D6FD}=1$ , 0.01 and 5, respectively (see figure 5 a,b). We also obtain the leading-order asymptotic behaviour for the particle velocity as

(4.8) $$\begin{eqnarray}\displaystyle q^{(2)}(0,t) & {\approx} & \displaystyle c(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FE})+\frac{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D6FC}^{2}(\unicode[STIX]{x03C0}t)^{1/2}}+O(t^{-3/2}),\quad t\gg 0,\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle c(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FE}) & = & \displaystyle \int _{-\infty }^{\infty }\frac{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FE}}{(\unicode[STIX]{x1D6FC}-k^{2})^{2}+k^{2}\unicode[STIX]{x1D6FE}^{2}}\,\text{d}k=\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D6FC}},\end{eqnarray}$$

where $c(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FE})$ is the velocity of the sedimenting particle at long times.

Figure 5.

(a) Particle velocity relative to fluid, $q^{(2)}(0,t)$ , along the $y^{(2)}$ -direction and (b) particle trajectory as a function of time for three density ratios, $\unicode[STIX]{x1D6FD}=0.01,1,5$ , corresponding to light, neutrally buoyant and heavy particle when acted upon by gravity for $S=0.01,~\unicode[STIX]{x1D70E}=-1$ in (4.7). The asymptotic expression for velocity in (4.9) is plotted as the dotted line in (a), and the dashed line in (b) is the trajectory of the particle without accounting for Basset history. The time taken by particles to reach $y^{(2)}=0$ when accounting for Basset history for these density ratios is $t=1.1,1.15,1.2$ and without Basset history it takes $t=1,1.01,1.04$ .

Note that the terminal velocity $c(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FE})$ is in fact the same quantity as that obtained for a particle in the absence of the Basset history force and, more importantly, this terminal velocity is independent of $\unicode[STIX]{x1D6FE}$ . The expression $\unicode[STIX]{x1D70E}\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}$ is the constant velocity obtained from a balance of gravitational forces and steady drag in our notation. The evaluation of the integral expression for $c(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FE})$ is via an appeal to the residue theorem, albeit for three different cases: $\unicode[STIX]{x1D6FE}^{2}/\unicode[STIX]{x1D6FC}<4$ , $\unicode[STIX]{x1D6FE}^{2}/\unicode[STIX]{x1D6FC}=4$ and $\unicode[STIX]{x1D6FE}^{2}/\unicode[STIX]{x1D6FC}>4$ . In all three cases, the integral has the same value.

We observe that a sedimenting particle attains its terminal velocity at a rate of $t^{-1/2}$ , which is far slower than the exponential trend as predicted without Basset history. It is worth noting that a sedimenting particle’s relaxation is even slower than the $t^{-3/2}$ trend observed for the force-free relaxation (4.3). Another interesting feature of this dynamics, which would be missed on neglecting the Basset history integral, is that the long-time trajectory of the particle is given by $y^{(2)}\sim -1/(y^{(1)}-\unicode[STIX]{x1D709})$ (with $\unicode[STIX]{x1D709}$ being the final horizontal coordinate). This behaviour is anomalous in the sense that it is neither projectile nor ballistic.

The previous remarks suggest that Basset history could play a significant role in complex flows. For example, the sedimentation rates of individual particles will determine how frequently particles collide with one another and coalesce as they sink. This question is of importance to estimating carbon sequestration by plankton in the ocean, for raindrop growth, and for pulverised coal in thermal power plants. To this end, here we make a simple estimate for how a particle grows in size as it sediments through a sea of smaller particles.

It is simpler in this case to work with the dimensional equation, obtained by replacing $\unicode[STIX]{x1D70E}\rightarrow (\unicode[STIX]{x1D6FD}-1)g/R$ and $S\rightarrow \unicode[STIX]{x1D70F}_{p}$ , where $g$ is gravity and $\unicode[STIX]{x1D70F}_{p}=a^{2}/\unicode[STIX]{x1D708}$ is the particle time scale, so (4.7) reads as

(4.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x03C0}q^{(2)}(0,t)=\int _{-\infty }^{\infty }\frac{\text{i}k\text{e}^{-k^{2}t}v_{0}^{(2)}}{\displaystyle \frac{1}{R\unicode[STIX]{x1D70F}_{p}}-k^{2}+\frac{\text{i}k}{R\sqrt{\unicode[STIX]{x1D70F}_{p}}}}\,\text{d}k-\int _{-\infty }^{\infty }\frac{\displaystyle \frac{(\unicode[STIX]{x1D6FD}-1)g}{R^{2}\sqrt{\unicode[STIX]{x1D70F}_{p}}}(1-\text{e}^{k^{2}t})}{\left(\displaystyle \frac{1}{R\unicode[STIX]{x1D70F}_{p}}-k^{2}\right)^{2}+\displaystyle \frac{k^{2}}{R^{2}\unicode[STIX]{x1D70F}_{p}}}\,\text{d}k. & & \displaystyle\end{eqnarray}$$

Now $k$ has the dimensions of $t^{-1/2}$ , and $v_{0}^{(2)}$ and $q^{(2)}(0,t)$ are dimensional velocities. We follow one large particle of growing radius $a$ as it sediments through a quiescent fluid with a uniform distribution of identical small particles of radius $a_{s}<a$ occupying a solid fraction $\unicode[STIX]{x1D711}$ . The larger particle falls faster than the smaller ones, and in a short time $\unicode[STIX]{x0394}t$ it collides with all small particles whose centres lie within a cylinder of radius $(a+a_{s})$ and height $q^{(2)}(0,t)\unicode[STIX]{x0394}t$ . Assuming that all small particles it collides with will coalesce with the bigger particle and that $a\gg a_{s}$ , its volume increases by $\unicode[STIX]{x1D711}\unicode[STIX]{x03C0}a^{2}q^{(2)}(0,t)\unicode[STIX]{x0394}t$ . Thus we get an evolution for the particle time scale as

(4.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{p}^{(n)}=\left(\unicode[STIX]{x1D70F}_{p}^{3/2}+\frac{3}{4}\unicode[STIX]{x1D711}\unicode[STIX]{x1D70F}_{p}\frac{q^{(2)}(0,t)\unicode[STIX]{x0394}t}{\sqrt{\unicode[STIX]{x1D708}}}\right)^{2/3}, & & \displaystyle\end{eqnarray}$$

which modifies $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FE}$ in (3.16). After each time interval $\unicode[STIX]{x0394}t$ , we reinitialise the system with a new initial condition $v_{0}^{(2)}$ , being the velocity from the previous time instant, and a new relaxation time $\unicode[STIX]{x1D70F}_{p}^{(n)}$ . For a droplet of $\unicode[STIX]{x1D708}\sim 10^{-5}~\text{m}^{2}~\text{s}^{-1}$ , $a\sim 0.3$  mm and $\unicode[STIX]{x1D6FD}=10$ , we get $\unicode[STIX]{x1D70F}_{p}\sim 10^{-2}$  s. In figure 6(a) we plot the dynamics of such a particle; and in figure 6(b) the particle size is represented in terms of its time scale $\unicode[STIX]{x1D70F}_{p}$ as a function of time for two scenarios: accounting for (i) history effects (solid lines) and (ii) only Stokes drag (dashed lines). We see that by neglecting the history force we will grossly overestimate particle size at a given time.

Figure 6.

(a) Velocity as a function of time of a particle growing in size as it sediments through a uniform distribution of smaller particles and coalesces with those it collides with. (b) The size of this particle as a function of time in terms of its relaxation time scale $\unicode[STIX]{x1D70F}_{p}$ with initial $\unicode[STIX]{x1D70F}_{p}=0.01$  s,  $v_{0}^{(2)}=0.01$ , $\unicode[STIX]{x1D708}=10^{-5}~\text{m}^{2}~\text{s}^{-1}$ and $\unicode[STIX]{x1D6FD}=10$ in an environment with solid fraction $\unicode[STIX]{x1D711}=0.002$ . Solid lines correspond to dynamics accounting for history effects, while dashed lines are for particle relaxation with Stokes drag alone.

We have considered only vertical trajectories in this example. However, in gravitational settling in turbulent particulate flows, the anomalous trajectory could introduce non-trivial changes in the dynamics, which will be a topic for future study.

4.3 Example 3: particle in an oscillatory background

We now consider the behaviour of a single particle with a background flow that oscillates with a single frequency. The forcing function is $f^{(2)}(s)=\sin (\unicode[STIX]{x1D706}s)$ and the first term in equation (3.16) can be rewritten as

(4.12) $$\begin{eqnarray}\displaystyle \mathscr{I} & = & \displaystyle \int _{-\infty }^{\infty }\int _{0}^{t}\frac{\text{i}k\text{e}^{-k^{2}m}f^{(2)}(t-m)}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\,\text{d}m\,\text{d}k,\quad \text{with}~m=(t-s),\end{eqnarray}$$

(4.13) $$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }\frac{\text{i}\unicode[STIX]{x1D706}k}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\left[\frac{\unicode[STIX]{x1D706}\text{e}^{-k^{2}t}-\unicode[STIX]{x1D706}\cos (\unicode[STIX]{x1D706}t)+k^{2}\sin (\unicode[STIX]{x1D706}t)}{\unicode[STIX]{x1D706}^{2}+k^{4}}\right]\text{d}k.\end{eqnarray}$$

The solution for the scenario with Stokes drag (but no Basset history integral) is

(4.14) $$\begin{eqnarray}\displaystyle q^{(2)}(0,t)=q^{(2)}(0,0)\text{e}^{-\unicode[STIX]{x1D6FC}t}+\frac{\unicode[STIX]{x1D706}\text{e}^{-\unicode[STIX]{x1D6FC}t}}{\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D706}^{2}}+\frac{\unicode[STIX]{x1D6FC}\sin (\unicode[STIX]{x1D706}t)-\unicode[STIX]{x1D706}\cos (\unicode[STIX]{x1D706}t)}{\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D706}^{2}}. & & \displaystyle\end{eqnarray}$$

We evaluate this velocity expression as a function of time in figure 7(a) and find that, with Basset history, particles oscillate with smaller amplitudes at early time, and that they attain their final periodic state faster than when Basset history is neglected. This behaviour, governed by the short-time dynamics, is qualitatively different from that of the cases we studied in earlier sections, where Basset history significantly slowed down the attainment of terminal velocity. In particular, the Basset history integral is not simply an additional drag. We see that expressions (4.13) and (4.14) asymptote to similar expressions, albeit with a phase difference. We can quantify the phase difference $\unicode[STIX]{x1D719}$ that persists at long times between a particle with history effects and with only Stokes drag at long times. This is shown in figure 7(b) for Stokes number in the range $S=0.1$ –1 and for three different density ratios, $\unicode[STIX]{x1D6FD}=0.01,1,5$ . We should however remind ourselves that, though we have chosen $\unicode[STIX]{x1D706}\sim O(1)$ where particle relaxation time scale and forcing are similar, the solution expression is general and applies to all scenarios with frequencies capped by the diffusive time scale.

Figure 7.

(a) Particle velocity relative to oscillating background flow when $\unicode[STIX]{x1D6FD}=5$ , $S=1$ for three different oscillating frequencies, $\unicode[STIX]{x1D706}=0.5,1,5$ . The dashed lines correspond to the solution with Stokes drag alone, and the solid line accounts for Basset history as well, from (4.13). (b) Phase difference $\unicode[STIX]{x1D719}$ between solutions including and excluding history effects as a function of Stokes number, $S$ , for a fixed $\unicode[STIX]{x1D706}=1$ and three density ratios, $\unicode[STIX]{x1D6FD}=0.01,1,5$ .

5.1 Example 4: migration in plane Couette flow

The Couette flow velocity profile is given by $u^{(1)}=\unicode[STIX]{x1D706}y^{(2)}$ . As before, the superscripts $(\cdot )^{(1)}$ and $(\cdot )^{(2)}$ represent the horizontal and vertical components, respectively. The particle velocity evolution in such a background flow is given by

(5.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\left[\begin{array}{@{}c@{}}q^{(1)}(0,t)\\ q^{(2)}(0,t)\end{array}\right]=\underbrace{\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D6FC} & \unicode[STIX]{x1D706}\\ 0 & \unicode[STIX]{x1D6FC}\end{array}\right]}_{\mathbb{M}}\left[\begin{array}{@{}c@{}}q^{(1)}(0,t)\\ q^{(2)}(0,t)\end{array}\right]+\unicode[STIX]{x1D6FE}\underbrace{\left[\begin{array}{@{}cc@{}}1 & 0\\ 0 & 1\end{array}\right]}_{\mathbb{I}}\left[\begin{array}{@{}c@{}}q_{x}^{(1)}(0,t)\\ q_{x}^{(2)}(0,t)\end{array}\right]. & & \displaystyle\end{eqnarray}$$

Since the boundary condition is linear with constant coefficients, we may follow a similar procedure as earlier. We now have two diffusion equations coupled at an interface, where the interfacial boundary condition in each is given by the MR equation (3.6). Note, however, that the diffusion equations themselves (i.e. in the bulk) are decoupled from each other. As a result, the global relation for this system of diffusion equations is just the vector version of (2.13):

(5.2) $$\begin{eqnarray}\displaystyle \text{e}^{\unicode[STIX]{x1D714}(k)T}\hat{\boldsymbol{q}}(k,T)+\boldsymbol{g}_{1}(k,T)+\text{i}k\boldsymbol{g}_{0}(k,T)=0,\quad k\in \mathbb{C}^{-}. & & \displaystyle\end{eqnarray}$$

The boundary condition becomes

(5.3) $$\begin{eqnarray}\displaystyle (\boldsymbol{q}(0,T)\text{e}^{\unicode[STIX]{x1D714}T}-\boldsymbol{q}(0,0))-k^{2}\boldsymbol{g}_{0}(k,T)=-\mathbb{M}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D706})\boldsymbol{g}_{0}(k,T)+\unicode[STIX]{x1D6FE}\boldsymbol{g}_{1}(k,T). & & \displaystyle\end{eqnarray}$$

Eliminating $\boldsymbol{g}_{1}(k,T)$ just as earlier from the above equation, we get for $\boldsymbol{g}_{0}(k,T)$

(5.4) $$\begin{eqnarray}\displaystyle ((k^{2}-\text{i}k\unicode[STIX]{x1D6FE})\mathbb{I}-\mathbb{M}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D706}))\boldsymbol{g}_{0}(k,T)=\{\boldsymbol{q}(0,T)\text{e}^{\unicode[STIX]{x1D714}T}-\boldsymbol{q}(0,0)+\unicode[STIX]{x1D6FE}\hat{\boldsymbol{q}}(k,T)\text{e}^{\unicode[STIX]{x1D714}T}\}, & & \displaystyle\end{eqnarray}$$

where

(5.5) $$\begin{eqnarray}\displaystyle ((k^{2}-\text{i}k\unicode[STIX]{x1D6FE})\mathbb{I}-\mathbb{M}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D706}))^{-1}=\left[\begin{array}{@{}cc@{}}(k^{2}-\unicode[STIX]{x1D6FC}-\text{i}k\unicode[STIX]{x1D6FE})^{-1} & \unicode[STIX]{x1D706}(k^{2}-\unicode[STIX]{x1D6FC}-\text{i}k\unicode[STIX]{x1D6FE})^{-2}\\ 0 & (k^{2}-\unicode[STIX]{x1D6FC}-\text{i}k\unicode[STIX]{x1D6FE})^{-1}\end{array}\right]. & & \displaystyle\end{eqnarray}$$

Multiplying the above equation by $\text{i}k\text{e}^{-k^{2}t}$ and integrating over $\unicode[STIX]{x2202}{\mathcal{D}}^{-}$ we obtain the solution

(5.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x03C0}\boldsymbol{q}(0,t)=-\int _{-\infty }^{\infty }\text{i}k\text{e}^{-k^{2}t}\left[\begin{array}{@{}cc@{}}(k^{2}-\unicode[STIX]{x1D6FC}-\text{i}k\unicode[STIX]{x1D6FE})^{-1} & \unicode[STIX]{x1D706}(k^{2}-\unicode[STIX]{x1D6FC}-\text{i}k\unicode[STIX]{x1D6FE})^{-2}\\ 0 & (k^{2}-\unicode[STIX]{x1D6FC}-\text{i}k\unicode[STIX]{x1D6FE})^{-1}\end{array}\right]\boldsymbol{q}(0,0)\,\text{d}k. & & \displaystyle\end{eqnarray}$$

On the other hand, the solution to the MR equation without the Basset history integral simplifies to

(5.7) $$\begin{eqnarray}\displaystyle \boldsymbol{q}(0,t)=\text{e}^{-\unicode[STIX]{x1D6FC}t}\left[\begin{array}{@{}cc@{}}1 & -\unicode[STIX]{x1D706}t\\ 0 & 1\end{array}\right]\boldsymbol{q}(0,0). & & \displaystyle\end{eqnarray}$$

The results are provided in figure 8. Figure 8(a) and (d), respectively, for Stokes number $S=0.01$ and 1, show particle trajectories with and without Basset history. The $y^{(2)}$ coordinate has saturated to the values shown. Clearly the migratory behaviour is different, and neglecting Basset history predicts the final locations of the particles incorrectly. Furthermore, the relaxation time is longer with Basset history, i.e. the velocity decay is slower, as seen in figure 8(b,e). Note that the relative velocity is plotted here, so each particle at long time attains the local velocity of the Couette flow ( $=\unicode[STIX]{x1D706}y^{(2)}$ ). The solution expressions (5.6)–(5.7) readily lead to the following expressions for the sensitivity of the solution to a particle’s initial velocity:

(5.8) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x1D6FF}\boldsymbol{q}(0,t)}{\unicode[STIX]{x1D6FF}\boldsymbol{q}(0,0)}\right)_{His}=-\frac{1}{\unicode[STIX]{x03C0}}\int _{-\infty }^{\infty }\text{i}k\text{e}^{-k^{2}t}\left[\begin{array}{@{}cc@{}}(k^{2}-\unicode[STIX]{x1D6FC}-\text{i}k\unicode[STIX]{x1D6FE})^{-1} & \unicode[STIX]{x1D706}(k^{2}-\unicode[STIX]{x1D6FC}-\text{i}k\unicode[STIX]{x1D6FE})^{-2}\\ 0 & (k^{2}-\unicode[STIX]{x1D6FC}-\text{i}k\unicode[STIX]{x1D6FE})^{-1}\end{array}\right]\,\text{d}k, & \displaystyle\end{eqnarray}$$

(5.9) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x1D6FF}\boldsymbol{q}(0,t)}{\unicode[STIX]{x1D6FF}\boldsymbol{q}(0,0)}\right)_{Sto}=\text{e}^{-\unicode[STIX]{x1D6FC}t}\left[\begin{array}{@{}cc@{}}1 & -\unicode[STIX]{x1D706}t\\ 0 & 1\end{array}\right]. & \displaystyle\end{eqnarray}$$

The eigenvalues of this sensitivity matrix provide the Lyapunov exponents, and both eigenvalues are equal to 1 for the scenarios with and without Basset history. However, the eigenvectors corresponding to these matrices are different and in figure 8(c,f) we plot the evolution of the off-diagonal component for two different $S=0.01,1$ .

Figure 8.

Particle dynamics in Couette flow. Solid lines are the solutions in (5.6), while dashed lines indicate dynamics upon neglecting Basset history in (5.7). (a–c) Corresponds to solution for $S=0.01$ and (d–f) for $S=1$ for a fixed density ratio of $\unicode[STIX]{x1D6FD}=5$ . (a,d) Horizontal component $q^{(1)}(0,t)$ and vertical component $q^{(2)}(0,t)$ , respectively, of particle velocity relative to fluid velocity with initial velocity of particle relative to fluid $\boldsymbol{q}(0,0)=1$ . (b,e) Particle trajectory as a function of time. We see that the particle velocity is the same in both directions, however decaying differently. (c,f) Off-diagonal component of sensitivity matrix in (5.9) which is due to shear in the flow.

5.2 Solution for nonlinear spatial flows

So far, we have been able to explicitly compute the particle velocity, from $\tilde{\boldsymbol{g}}_{0}(\unicode[STIX]{x1D714},T)$ ; see (3.14). However, when the spatial dependence of the fluid flow is nonlinear, i.e. when $\boldsymbol{f}$ depends on $\boldsymbol{y}$ , equation (3.14) is not an explicit expression for the particle velocity. Consequently we resort to a numerical scheme to evolve the dynamics of the particle, which we present in this section. Our numerical method has the following advantages:

1. (i) The method completely eliminates the need to store the entire trajectory of the particle to evaluate the Basset history term.

2. (ii) Our method achieves spectral accuracy by employing Chebyshev polynomials to evaluate the integrals involved.

3. (iii) The method involves a fixed memory cost (of $O(N)$ , $N$ being the number of Chebyshev modes) independent of time, making long-time simulations no more expensive than short-time ones. This is useful to compute statistics in turbulent simulations.

4. (iv) The effect of history can be recorded (with fixed memory size), enabling simulations to be restarted. Such a restart is impossible in other techniques owing to exponentially large memory costs in evaluating the integral. While preparing the paper we became aware of the work of Parmar et al. (Reference Parmar, Annamalai, Balachandar and Prosperetti2018) who approximate the history kernel by a summation of exponentials and achieve the desired accuracy. However, we would like to highlight that our technique does not make any approximation to the kernel and we show later that all information due to history can be stored in a variable whose computation is of fixed memory throughout particle evolution. Moreover, the form of the solution in (3.16) is specially useful for studying asymptotic behaviour (authors’ work in preparation), unlike other techniques.

In the following subsection we derive the equations we ultimately solve for the particle trajectory and in the subsequent subsection we provide the details of the numerical implementation.

5.2.1 Evolution equation for Basset history integral

Figure 9.

Schematic showing numerical scheme where we evaluate $\hat{\boldsymbol{q}}(k,t)$ at intervals $t_{0},~(t_{0}+\unicode[STIX]{x0394}t)$ from the solution at $t\in [t_{0}-\unicode[STIX]{x0394}t,t_{0}],~[t_{0},t_{0}+\unicode[STIX]{x0394}t]$ and so on in order to eliminate the need to store $\boldsymbol{q}(0,t)$ from $t\in [0,t_{0}+\unicode[STIX]{x0394}t]$ .

We begin by dividing the $x$ – $t$ plane into domains of size $[0,\infty )\times [t_{0}-\unicode[STIX]{x0394}t,t_{0}]$ , where $t_{0}=m\unicode[STIX]{x0394}t$ for some integer $m$ , as shown in figure 9. Our procedure involves two steps, repeated for each such domain:

1. (1) Using the initial condition $\boldsymbol{q}(x,t_{0}-\unicode[STIX]{x0394}t)$ in each domain (bottom wall of that domain in figure 9), we compute the boundary condition $\boldsymbol{q}(0,t),~t\in (t_{0}-\unicode[STIX]{x0394}t,t_{0}]$ (left wall of that domain). This is done by solving the MR equation in the form given by (5.11) below. This equation is the extension of (3.16) to the case of arbitrary initial conditions. Moreover, the forcing $\boldsymbol{f}$ is now a function of $\boldsymbol{q}(0,t)$ and/or $\boldsymbol{y}(t)$ . So we start with a guess for $\boldsymbol{q}(0,t)$ and $\boldsymbol{y}(t)$ , and iterate using a nonlinear solver. We thus have the Dirichlet condition on the left wall of figure 9, over $[t_{0}-\unicode[STIX]{x0394}t,t_{0}]$ .

2. (2) Using the computed boundary condition $\boldsymbol{q}(0,t)$ from the previous step and the known initial condition $\boldsymbol{q}(x,t_{0}-\unicode[STIX]{x0394}t)$ , we directly compute $\boldsymbol{q}(x,t_{0})$ , the top wall of that domain, by using (5.13), which will be derived below. This step is equivalent to solving the diffusion equation with a Dirichlet boundary condition, and the solution of this boundary-value problem leads to the initial condition for the next domain, $t\in (t_{0},t_{0}+\unicode[STIX]{x0394}t]$ . Through an efficient use of the global relation, we do not require the explicit solution expression $\boldsymbol{q}(x,t_{0})$ , but only the Fourier transform $\hat{\boldsymbol{q}}(k,t_{0})$ , which simplifies the calculations considerably.

We recognise that our system size (namely, the number of dependent variables) is higher for a nonlinear case, and this is the price we need to pay to eliminate the rising memory costs coming from the non-local time integral term.

In the first domain alone $[0,\infty )\times [0,\unicode[STIX]{x0394}t]$ , equation (3.16) may be used to obtain the particle velocity. We recognise that $\boldsymbol{f}$ now depends on $\boldsymbol{y}$ and $\boldsymbol{q}(0,t)$ . Along with (3.7), where we set $T=\unicode[STIX]{x0394}t$ , we now have two equations for two unknowns, $\boldsymbol{y}$ and $\boldsymbol{q}(0,t)$ , over the time interval $[0,\unicode[STIX]{x0394}t]$ . In a standard numerical method, one would approximate the first integral of (3.16) using a quadrature rule. We, however, adopt an alternative route (partly to ensure accuracy) by representing the unknowns $\boldsymbol{q}(0,t)$ and $\boldsymbol{y}(t)$ in terms of Chebyshev polynomials over the interval $[0,\unicode[STIX]{x0394}t]$ . All integrals are readily evaluated to a high degree of accuracy. With the Chebyshev polynomial representation, equations (3.16) and (3.7) are solved using a nonlinear solver, here Newton’s method, for the coefficients of the Chebyshev polynomials. Once we have $\boldsymbol{q}(0,t)$ (in terms of Chebyshev polynomials), we solve the diffusion equation with this Dirichlet condition to obtain the solution $\boldsymbol{q}(x,\unicode[STIX]{x0394}t)$ . The evaluation of such new initial conditions at every $(t_{0}=m\unicode[STIX]{x0394}t)$ is possible precisely because of the relationship between the Basset history integral and the Neumann condition to the diffusion equation.

In all later domains, the initial condition $\boldsymbol{q}(x,t_{0})$ is no longer zero. Indeed $\boldsymbol{q}(x,t_{0})$ is the solution of the diffusion equation at this time. The global relation of the diffusion equation in the time interval $t\in (t_{0},t_{0}+\unicode[STIX]{x0394}t]$ can be written as

(5.10) $$\begin{eqnarray}\displaystyle & & \displaystyle \text{e}^{k^{2}(t_{0}+\unicode[STIX]{x0394}t)}\hat{\boldsymbol{q}}(k,(t_{0}+\unicode[STIX]{x0394}t))\nonumber\\ \displaystyle & & \displaystyle \quad =\text{e}^{k^{2}t_{0}}\hat{\boldsymbol{q}}(k,t_{0})-\text{i}k\tilde{\boldsymbol{g}}_{0}(k^{2},t_{0},t_{0}+\unicode[STIX]{x0394}t)-\tilde{\boldsymbol{g}}_{1}(k^{2},t_{0},t_{0}+\unicode[STIX]{x0394}t),\quad k\in \mathbb{C}^{-}.\end{eqnarray}$$

For the boundary terms represented by $\tilde{\boldsymbol{g}}_{i}(k^{2},t_{1},t_{2})$ , we have introduced the additional arguments $t_{1}$ and $t_{2}$ , which are the limits of the integral defining the time transform. This notation is also used in the forcing, as $\tilde{\boldsymbol{f}}(k^{2},t_{1},t_{2})$ . Following a procedure similar to that employed to derive (3.14), we obtain the following expression for $\boldsymbol{q}(0,t)$

(5.11) $$\begin{eqnarray}\displaystyle -\frac{\unicode[STIX]{x03C0}}{2}\boldsymbol{q}(0,t) & = & \displaystyle \underbrace{\int _{0}^{\infty }k\text{e}^{-k^{2}(t-t_{0})}\text{Im}\{\mathscr{H}(k,t_{0})\}\,\text{d}k}_{{\mathcal{H}}(t)}\nonumber\\ \displaystyle & & \displaystyle +\,\underbrace{\int _{0}^{\infty }\text{Im}\left\{\frac{k\text{e}^{-k^{2}t}\tilde{\boldsymbol{f}}(k^{2},t_{0},t)}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\right\}\,\text{d}k}_{{\mathcal{F}}(t)},\quad t\in (t_{0},t_{0}+\unicode[STIX]{x0394}t].\end{eqnarray}$$

The term $\mathscr{H}(k,t_{0})$ in the integrand of the first integral in the expression above is the modified initial condition given by

(5.12) $$\begin{eqnarray}\displaystyle \mathscr{H}(k,t_{0})=\left\{\frac{\boldsymbol{q}(0,t_{0})+\unicode[STIX]{x1D6FE}\hat{\boldsymbol{q}}(k,t_{0})}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}\right\}. & & \displaystyle\end{eqnarray}$$

This term encodes all the information about the history of the particle from $t\in (0,t_{0}]$ . Equation (5.11) is thus the extension of (3.16) for a general initial condition to the diffusion equation. The renormalised initial condition $\mathscr{H}(k,t)$ itself evolves in time and this evolution is once again obtained from the global relation; see appendix B for details. We state the relevant equation here:

(5.13) $$\begin{eqnarray}\displaystyle \mathscr{H}(k,t_{0})=\text{e}^{-k^{2}\unicode[STIX]{x0394}t}\mathscr{H}(k,t_{0}-\unicode[STIX]{x0394}t)-\underbrace{\text{e}^{-k^{2}t_{0}}\int _{t_{0}-\unicode[STIX]{x0394}t}^{t_{0}}\text{e}^{k^{2}s}\boldsymbol{q}(0,s)\,\text{d}s}_{\mathscr{I}_{1}}-\underbrace{\text{e}^{-k^{2}t_{0}}\frac{\,\tilde{\boldsymbol{f}}(k^{2},t_{0}-\unicode[STIX]{x0394}t,t_{0})}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}}}_{\mathscr{I}_{2}}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Note that all terms on the right-hand side of (5.13), including the forcing term, are known since they depend on the previous time interval $[t_{0}-\unicode[STIX]{x0394}t,t_{0}]$ . Equation (5.13) may also be written as a dynamical system for individual modes of history as

(5.14) $$\begin{eqnarray}\displaystyle \dot{\mathscr{H}}(k,t)+k^{2}\mathscr{H}(k,t)=\boldsymbol{q}(0,t)+\frac{\,\boldsymbol{f}(t,\boldsymbol{y}(t),\boldsymbol{q}(0,t))}{\unicode[STIX]{x1D6FC}-k^{2}+\text{i}k\unicode[STIX]{x1D6FE}},\quad t\in [t_{0}-\unicode[STIX]{x0394}t,t_{0}]. & & \displaystyle\end{eqnarray}$$

Equations (5.11) and (5.13) represent the solution to the full MR equations for a particle in a time-dependent inhomogeneous flow field. Note that the individual modes $k$ of $\mathscr{H}(k,t)$ in (5.13) evolve independently. However, the expression for the velocity in (5.11) depends on all of them, and indeed depends on the integral of the history.

We are now enabled to march repeatedly till any $t_{0}$ from any $(t_{0}-\unicode[STIX]{x0394}t)$ using the set of initial conditions $\hat{\boldsymbol{q}}(k,t_{0}-\unicode[STIX]{x0394}t)$ , solving equation (5.11) for $\boldsymbol{y}(t)$ and $\boldsymbol{q}(0,t)$ in the interval $(t_{0}-\unicode[STIX]{x0394}t)$ to $t_{0}$ , and using these Dirichlet boundary conditions in (5.13) to get $\hat{\boldsymbol{q}}(k,t_{0})$ . Once again, we emphasise that the above procedure eliminates the need to store the entire particle velocity from the initial time instant. The effect of the Basset history integral is encoded in $\mathscr{H}(k,t)$ .

5.3 Numerical technique

We now detail a numerical scheme to solve equations (5.11)–(5.13). The continuous variables $t$ and $k$ are restricted to a set of discrete values over their respective intervals, $t\in (t_{0}-\unicode[STIX]{x0394}t,t_{0}]$ and $k\in [0,\infty )$ , leading to discrete times $t_{j}$ and discrete wavenumbers $k_{l}$ where $j=1,\ldots ,L$ and $l=0,\ldots ,N$ . In the case of $t$ we employ Chebyshev nodes under an affine transform, whereas for $k$ we employ a mapped Chebyshev node appropriate for rational Chebyshev approximations. We then expand $\boldsymbol{q}(0,t)$ and $\mathscr{H}(k,t_{0})$ in a Chebyshev polynomial basis $T_{n}(t)$ (rational Chebyshev polynomials for $\mathscr{H}$ ). We adopt a collocation approach to solve the resulting nonlinear equation.

The discrete version of (5.11) can be written as a nonlinear equation for the coefficients as

(5.15) $$\begin{eqnarray}\displaystyle \mathscr{T}(t_{j})\equiv \frac{\unicode[STIX]{x03C0}}{2}\boldsymbol{q}(0,t_{j})+{\mathcal{H}}(t_{j})+{\mathcal{F}}(t_{j})=0,\quad t_{j}=t_{0}+\frac{\unicode[STIX]{x0394}t}{2}\left(1+\cos \left(\frac{j\unicode[STIX]{x03C0}}{L}\right)\right), & & \displaystyle\end{eqnarray}$$

where $j=1,\ldots ,L$ . Details on how to compute each term in the above equation are presented in appendix C. The contribution from $\mathscr{H}(k,t_{0})$ is captured in ${\mathcal{H}}(t_{j})$ , whereas the contribution from the nonlinear forcing due to the fluid velocity field is present in ${\mathcal{F}}(t_{j})$ . An equation for the position $\boldsymbol{y}$ can be obtained from the discretised version (3.7). A standard Newton root finder is employed to solve for $\boldsymbol{q}(0,t)$ and $\boldsymbol{y}(t)$ in the time interval $(n\unicode[STIX]{x0394}t,(n+1)\unicode[STIX]{x0394}t]$ for each $n=0,1,\ldots$ with an appropriate initial guess, typically $\boldsymbol{q}(0,t_{j})=-2/\unicode[STIX]{x03C0}{\mathcal{H}}(t_{j})$ . Equation (5.13) on the other hand is a straightforward integral evaluation using known values of $\boldsymbol{q}(0,t_{j})$ and $\mathscr{H}(k_{j},t_{0}-\unicode[STIX]{x0394}t)$ from the previous time step. We show in the appendix that expanding $\boldsymbol{q}(0,t)$ and $\mathscr{H}(k,t_{0}-\unicode[STIX]{x0394}t)$ in Chebyshev basis allows us to precompute several terms, which speeds up evaluation of expressions in (5.13).

Figure 10.

Error estimates for the solution to a relaxing particle using the numerical scheme developed in § 5.3. Error as a function of time for (a) different values of $N$ at fixed $L=20$ and (c) different values of $L$ at fixed $N=50,70$ . (b,d) Maximum of the error as we increase the number of Chebyshev modes in $k$ and  $t$ .

5.4 Convergence properties and memory requirement of the numerical scheme

To study convergence of our numerical scheme, we compare the exact velocity of a relaxing particle (for which we have an analytic expression) to the numerical solution obtained from the scheme described above. We first check the numerical method’s accuracy as a function of the number of Chebyshev nodes $N$ in wavenumber $\tilde{k}$ , where $\tilde{k}:~k\in [0,\infty ]\rightarrow [-1,1]$ . We also check the accuracy for number of Chebyshev nodes $L$ in $t\in [t_{0},t_{0}+\unicode[STIX]{x0394}t]$ within a time interval of length $\unicode[STIX]{x0394}t$ . The error in the calculated solution is given by

(5.16) $$\begin{eqnarray}\displaystyle \mathscr{E}(t_{j})=|q^{(1)}(0,t_{j})-q_{e}^{(1)}(0,t_{j})|, & & \displaystyle\end{eqnarray}$$

where $q_{e}^{(1)}(0,t_{j})$ is the exact solution evaluated at discrete $t_{j}$ and $q^{(1)}(0,t_{j})$ is the solution evaluated using the present numerical technique. In figure 10(a) we plot the error as a function of time for all the number of $\tilde{k}$ chosen and we see its magnitude decrease with increase in $N$ . In figure 10(b) we plot the maximum value of $\mathscr{E}(t_{j})$ . We clearly see that the error is small and decreases exponentially with increase in $N$ for fixed $L=20$ . We can reach a desired accuracy by choosing the appropriate number of Chebyshev nodes $N$ . Moving now to the convergence of the solution with increase in nodes $L$ in $\unicode[STIX]{x0394}t$ , we see there is a general increase in error in figure 10(d) for both fixed $N=50$ and 70. Such a behaviour is expected and not an anomaly. Firstly observe that the maximum error is attained, for all different values of $L$ and $N$ , at the first grid point, i.e.  $t_{1}$ in $t_{j}\in (0,\unicode[STIX]{x0394}t]$ . This is because of the $t^{-1/2}$ singularity at $t=0$ . As we increase the number of Chebyshev nodes in the interval of length $\unicode[STIX]{x0394}t$ , the first Chebyshev grid point moves closer to this singularity. We can easily account for this behaviour by increasing the number of modes in space $N$ . We show this by increasing the spatial Chebyshev modes from $N=50$ to 70 in figure 10(d), which results in a drop in the maximum error.

We emphasise the advantage that a Chebyshev basis provides in terms of computational efficiency. In order to compute the Chebyshev coefficients via the discrete Chebyshev transform (DCT), we only require $O(N\log N)$ operations, where $N$ is the number of Chebyshev modes. This is done by leveraging a fast Fourier transform (FFT) package to evaluate the DCTs involved in the process of computing Chebyshev coefficients.

The memory requirement for our technique comes from storing the variables $\boldsymbol{q}(0,t_{j})$ and $\mathscr{H}(k_{j},t_{0})$ . Considering $L$ modes in $t_{j}$ , $N$ modes in $\tilde{k}_{j}$ demands only storage of $O(N+L)$ values for all times, whereas the other techniques require a linearly increasing memory cost of $O(t/\unicode[STIX]{x0394}t)$ . For $N=30,~L=20$ , one needs to store only 50 values, which remain fixed throughout the simulation. In this estimate, of course, we do not account for the precomputed matrices, which turn out to be three $N\times N$ matrices (see appendix C) resulting in another 2700 values for $N=30$ , but this initial cost remains constant and helps in shorter computation time.

5.5 Example 5: single point vortex and caustics

The stationary flow field corresponding to a single vortex is

(5.17) $$\begin{eqnarray}\displaystyle \boldsymbol{u}=\frac{\unicode[STIX]{x1D6E4}}{2\unicode[STIX]{x03C0}}\hat{\boldsymbol{e}}_{z}\times \frac{\,\boldsymbol{y}(t)}{|\,\boldsymbol{y}|^{2}}. & & \displaystyle\end{eqnarray}$$

Though this flow field does not have a time-varying component, it still has contributions from both $(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}$ and $(\boldsymbol{q}(0,t)\cdot \unicode[STIX]{x1D735})\boldsymbol{u}$ . This makes the system of evolution equations for particle position $\boldsymbol{y}(t)$ and velocity $\boldsymbol{q}(0,t)$ coupled and nonlinear. The solution is given by

(5.18) $$\begin{eqnarray}\displaystyle \boldsymbol{y}(t)=\boldsymbol{y}(0)+\int _{0}^{t}\boldsymbol{q}(0,s)\,\text{d}s+\int _{0}^{t}\frac{\unicode[STIX]{x1D6E4}}{2\unicode[STIX]{x03C0}}\hat{\boldsymbol{e}}_{z}\times \frac{\,\boldsymbol{y}(s)}{|\,\boldsymbol{y}|^{2}}\,\text{d}s, & & \displaystyle\end{eqnarray}$$

for the velocity evolution

(5.19) $$\begin{eqnarray}\displaystyle \boldsymbol{q}_{t}(0,t)+\unicode[STIX]{x1D6FC}\boldsymbol{q}(0,t)-\unicode[STIX]{x1D6FE}\boldsymbol{q}_{x}(0,t)=\left(\frac{1}{R}-1\right)(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}-(\boldsymbol{q}(0,t)\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}. & & \displaystyle\end{eqnarray}$$

Using the numerical technique described in the previous section, we solve for position and velocity for a Stokes number of $S=0.01$ and density ratio $\unicode[STIX]{x1D6FD}=20$ , and plot the radial and angular velocity in figure 11(a,b). Different lines in the figure correspond to different initial position separated by a distance of 0.15 though all particles start with an initial velocity of $(0.05,0.05)$ . Once the particles get centrifuged far away from the vortex, the trajectory follows $y^{(r)}(t)\sim t^{1/4}$ , where $y(r)$ is the radial distance from the point vortex. Thus the radial velocity $q^{(r)}(0,t)\sim t^{-3/4}$ , as in figure 11(c). This long-time behaviour is the same as that obtained by Ravichandran & Govindarajan (Reference Ravichandran and Govindarajan2015) without Basset history. The history integral affects the transient behaviour and the approach to the long-time limit.

Figure 11.

(a) Radial velocity, (b) angular velocity and (c) radial distance travelled by particles starting at five different locations separated by a non-dimensional distance of 0.15 at $S=0.01$ and $\unicode[STIX]{x1D6FD}=10$ , with an initial velocity of $\boldsymbol{v}_{0}=(0.05,0.05)$ and $\unicode[STIX]{x1D6E4}=2\unicode[STIX]{x03C0}$ . The solid black lines correspond to power laws of $t^{-3/4}$ or $t^{1/4}$ , as marked. Solid coloured lines represent simulations including the Basset history integral, Stokes drag and added-mass effects, while the dashed coloured lines are for Stokes drag with added mass.

We discussed in the introduction about caustics near a vortex. As shown analytically by Ravichandran & Govindarajan (Reference Ravichandran and Govindarajan2015), and in 2-D simulations of turbulence, particles inside a critical radius $r_{c}$ are evacuated rapidly, resulting in a region of high particle concentration in the neighbourhood of $r_{c}$ . In order to identify this critical radius, including the effects of Basset history, we track the radial separation $\unicode[STIX]{x0394}r$ between two particles that are initially placed close to each other. This is shown in figure 12(a) for various initial radial distances $r_{i}$ from the vortex of the inner particle (the one closer to the vortex). When $\unicode[STIX]{x0394}r$ crosses zero, caustics have occurred. We find that caustics occur for $r_{i}\leqslant 0.44$ , i.e.  $r_{c}=0.44$ , for this Stokes number $S=0.02$ . Also, the trajectories of particles that exhibit caustics take different paths compared to ones that do not, in the phase plane defined by $\unicode[STIX]{x0394}v_{r}$ versus $\unicode[STIX]{x0394}r$ , where $\unicode[STIX]{x0394}v_{r}$ is the difference in radial velocities for two closely placed particles (see Ravichandran & Govindarajan Reference Ravichandran and Govindarajan2015). This is shown in figure 12(b), where a negative value of $\unicode[STIX]{x0394}v_{r}$ implies that the inner particle travels faster than its neighbour. When $\unicode[STIX]{x0394}r$ becomes negative it has overtaken its neighbour. Thus trajectories that cross the $y$ -axis of the phase plane form caustics.

In figure 12(c) we plot the critical radius $r_{c}$ (found from the location where $\min (\unicode[STIX]{x0394}r(t))=0$ ) as a function of the Stokes number $S$ . Simulations with history integral and added-mass effects included are shown as solid lines and with just Stokes drag as dotted lines. The history integral and added mass together reduce the region over which caustics are formed, and may cause a reduction in the number densities of particles predicted without the history term. The implications for particle dynamics in turbulent fluid flows are beyond the realm of the present article and will be explored in our upcoming study.

Figure 12.

(a) Evolution of the difference in radial distance between particles $\unicode[STIX]{x0394}r$ for different initial particle locations $r_{i}$ separated by a distance $0.05$ for Stokes number $S=0.02$ . The critical particle location at which $\min (\unicode[STIX]{x0394}r)$ crosses the $t$ -axis is the location of the end of caustic formation. (b) The phase portrait defined by $\unicode[STIX]{x0394}v_{r}$ , the difference in radial velocity versus $\unicode[STIX]{x0394}r$ , also shows the transition from caustic formation to no caustic formation with increase in $r_{i}$ . (c) Critical radius $r_{c}$ for caustic formation as a function of Stokes number $S$ with history integral, added mass (solid line) and with Stokes drag alone (dashed line).