\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx,subfig}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2015 (2015), No. 104, pp. 1--19.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2015 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2015/104\hfil Three models for rectilinear particle motion]
{Three models for rectilinear particle motion with the Basset
history force}

\author[S. Xu, A. Nadim \hfil EJDE-2015/104\hfilneg]
{Shujing Xu, Ali Nadim}

\address{Shujing Xu \newline
Institute of Mathematical Sciences,
Claremont Graduate University,
Claremont, CA 91711, USA}
\email{flora.xushujing@gmail.com}

\address{Ali Nadim \newline
Institute of Mathematical Sciences,
Claremont Graduate University,
Claremont, CA 91711, USA}
\email{ali.nadim@cgu.edu}


\thanks{Submitted November 7, 2014. Published April 21, 2015.}
\subjclass[2010]{76T20, 45J05, 65L05}
\keywords{Particle motion; history force; integro-differential equations;
\hfill\break\indent Laplace transforms; numerical methods; multiphase flow}

\begin{abstract}
 We consider three model problems that describe rectilinear particle motion
 in a viscous fluid under the influence of the Basset history force.
 These problems consist of sedimentation starting from rest, impulsive
 motion in a quiescent fluid, and oscillatory sliding motion.
 The equations of motion are integro-differential equations with a weakly
 singular kernel. We derive analytical solutions to all three problems
 using Laplace transforms and discuss the mathematical relation between the
 sedimentation and impulsive start problems. We also compare several numerical
 schemes for solving the integro-differential equations and benchmark them
 against the analytical results.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\allowdisplaybreaks


\section{Introduction} \label{sec:Introduction}

The dynamics of a solid particle moving in a flowing liquid or gas is of great
interest in the field of fluid dynamics for it has broad application in physics,
biology, multiphase flow, etc. In the present paper, three fundamental
one-dimensional motions are examined in the presence of the history force:
sedimentation of a particle in a quiescent fluid in which the particle is
released from rest, impulsive motion where the particle is given an initial
velocity or impulse and subsequently relaxes under the influence of viscosity,
and an oscillatory sliding motion where a fluid-filled cartridge containing
the particle undergoes a back and forth motion.

For a small, spherical, non-deformable particle moving in an unbounded flow domain,
the equation of motion is given by \cite{maxey:883}
\begin{equation}
\begin{aligned}
m_p \frac{d\mathbf{V}}{dt}
&=m_f \frac{D\mathbf{u}}{Dt}
- \frac{1}{2} m_f \Bigl(\frac{d\mathbf{V}}{dt} - \frac{D\mathbf{u}}{Dt}\Bigr)
- 6\pi \mu a (\mathbf{V} - \mathbf{u})\\
&\quad + (m_p - m_f)\mathbf{g}
 - 6a^2 \sqrt{\pi \mu \rho_f} \int_{-\infty}^{t}\frac{1}{\sqrt{t - \tau}}
\Bigl(\frac{d  \mathbf{V}}{d\tau} - \frac{D\mathbf{u}}{D \tau} \Bigr) d\tau\,,
\end{aligned}\label{eqn:eom}
\end{equation}
where
$m_p$ is the particle mass;
$m_f$  the displaced fluid mass;
$\mathbf{V}(t)$  the particle velocity;
$\mathbf{u}(\mathbf{x}, t)$  the fluid velocity field;
$a$  the radius of the particle;
$\mu$  the dynamic viscosity of the fluid,
and $\rho_f$ its density.
Two time derivatives, $D/Dt$ and $d/dt$, are used to measure the rate-of-change
of fluid velocity seen by an observer moving along with the fluid and the particle,
respectively.

In addition to the pressure stress, added mass, viscous drag, and gravity as shown
on the right-hand side of \eqref{eqn:eom},  a history-type force depending
on the entire history of the motion is also known to play a role (the last term
in the equation). This Boussinesq-Basset force accounts for the temporal delay
in boundary layer development with changing relative velocity of bodies moving
through a fluid. This unsteady force acting on a particle submerged in a fluid
was first investigated by Joseph Valentin Boussinesq \cite{Boussinesq1885}
and Alfred Barnard Basset \cite{MR0134559}. It is also known simply as the
history force, a combination of both viscous and inertial contributions to the
force in that it depends on both the viscosity and density of the fluid and the
acceleration of the particle \cite{LovBrady1993}.
For the motion of rigid spheres in a viscous fluid, the regimes where the history
force might be negligible or is significant have been investigated. Parameters
that determine the significance of the history term include velocity fluctuations
 in time, contrasting fluid-to-particle density ratio, and the particle Reynolds
number \cite{  
FLM:353169, FLM:340128,  FLM:14341, ling:113301, 
LovBrady1993, lukerchenko2010basset, CambridgeJournals:397249, 
citeulike:11142061, reeks:1573,  visitskii2006,
Vojir1994547}.

The history term turns the equation of motion into an integro-differential
equation implicit in the dependent variable $\mathbf{V}(t)$ for it also appears
 inside the integrand and makes it challenging to solve for the motion.
As its name suggests, the history force involves an integral based on the
entire history of the motion up to the present time. For practical reasons,
the starting point in time is usually set as 0 instead of $-\infty$, under
the assumption that there is no prehistory before time 0. Even with such
simplification, for each time step, it is still necessary to make use of all
former values, which increases both the computational time and the memory
 requirements.  Another challenge lies in the fact that the integrand has a
singularity near its upper limit, which prevents one from using traditional
methods for performing the numerical quadrature. Overcoming these numerical
challenges is one aspect of the present work.

For certain problems involving the history force, exact analytical solutions
can be obtained. For these, Laplace transforms turn out to be the most fruitful
approach for obtaining the solutions in the time domain: e.g., for a sphere
undergoing free fall in a still, viscous fluid \cite{Brush1964}, a sphere in
creeping flow \cite{michaelides:1579}, and an isolated denser particle dropped
in the core of a vertical vortex \cite{candelier:1765}. 
Lovalenti and Brady  \cite{LovBrady1993}
pointed out one of the limitations of using Laplace transforms (in time) when
dealing with time-dependent coefficients. For a particle suspended in homogeneous
turbulence responding to the random velocity field, 
Mei et al  \cite{CambridgeJournals:397249}
resorted to Fourier transforms to obtain the analytical solution in the frequency
domain. \cite{yang:1822} evaluated the trajectories of an accelerating spherical
drop at low Reynolds number by transforming the result from the Fourier-transform
domain. Other analytical approaches include the use of Abel's theorem to
investigate linear models for a sphere falling through a Newtonian fluid 
\cite{Belmonte2000}; manipulating the operators, for example, writing the
integral
\[
I(f(t)) = \frac{1}{\sqrt\pi} \int_0^t \frac{f(x)}{\sqrt{t - x}} dx
= \frac{2}{\sqrt{\pi}} \int_0^{\sqrt t} f(t - x^2) dx
\]
 when the Cauchy problem
 is involved \cite{Visitskii2009548, Vodopaoyanov2010}, and the use of
the $\Gamma$ function similarly in \cite{yih1977fluid}. Coimbra et al 
\cite{MR1647206} generalized how to obtain the solution of 
the particle equation analytically for unsteady Stokes flows. 
First, they applied a fractional-differential
operator to the first-order integro-differential equation of motion in order
to transform the original equation into a second-order non-homogeneous equation,
and then they solved this last equation by the method of variation of
parameters.

On the numerical side, Brush et al \cite{Brush1964} proposed a trapezoidal-based method
under the assumption that on a small enough time interval, the derivatives in
the integrand could be roughly regarded as a constant and be brought outside
the integral sign, with the rest of the integrand evaluated analytically.
Taylor expansions were employed in \cite{Daitche2012,Daitche2011} to improve
the accuracy; however, this came at a cost of much longer computational
time. Bombardelli et al \cite{Bombardelli2008} considered the integral as a fractional
derivative that could be approximated by a summation series which had a first
order temporal accuracy. An approximation based on exponential functions was
described in \cite{hinsberg} and shown to have second-order accuracy; however,
the choice of exponential functions varies from one situation to the next.

In this work, we focus on three model problems whose exact solutions can be
obtained analytically. Among these, the exact solution for particle trajectories
in an oscillatory sliding motion is new. We also use these exact solutions to
benchmark some numerical schemes for solving such problems. The rest of the
present paper is structured as follows.
Section~\ref{sec:NumApp} introduces two numerical treatments for the
problematic history integral in detail. This is followed in
Section~\ref{sec:PhyApp} with the three physical problems involving
one-dimensional motion: sedimentation, impulsive motion, and oscillatory
sliding motion. For each problem, the scaled equation of motion is presented,
along with its analytical solution which is used to benchmark the corresponding
numerical methods. Selected numerical methods from the existing literature are
also applied to the same problems by way of comparison. Section~\ref{sec:con}
provides our summary and conclusions.

\section{Numerical approximation of the singular integrals}
\label{sec:NumApp}

For both methods that will be presented, the uniform time step is denoted
by $\Delta t$ in the temporal discretization and the initial time $t_1$
is taken to be 0. Thus, any given time $t_n$ is given by $t_n = (n-1)\Delta t$,
and the discrete analog of any function $f(t)$ at $t_n$ is denoted by
$f(t_n)$ or simply $f_n$.

\subsection{Forward difference hybrid}
\label{subsec: m1}
The first method introduced is a combination of trapezoidal rule and integration
 by parts carried out in two stages. It can be viewed as a refined version
of the trapezoidal-based method proposed by Brush in which it was assumed
that the integrand on a small interval is a constant. To introduce this two-stage
scheme, let us start with rewriting the integral history term in the form
\begin{equation}
I(t) = \int_0^t  \frac{f(\tau)}{\sqrt{t - \tau}} \,  d\tau \,.
\label{intHisVer01}
\end{equation}
In \eqref{eqn:eom}, $f \equiv  {d\mathbf{V}}/{d\tau} - {D\mathbf{u}}/{D\tau}$,
and is differentiable everywhere. When $\tau$ gets close to the upper limit $t$,
the singularity prevents us from applying the trapezoidal rule directly.
An intuitive idea is to isolate the troublemaker, i.e., break up the integral
into two parts, one containing the singularity on a very small interval
(as $\Delta t \to 0$) and the other being free of singularity.
That is, at any given time $t_n$, the integral can be evaluated by
$$
 I(t_n)=\int_0^{t_n} \frac{f(\tau)}{\sqrt{t_n - \tau}} d\tau
 = \int_0^{t_{n-1}} \frac{f(\tau)}{\sqrt{t_n - \tau}} d\tau
 + \int_{t_{n-1}}^{t_n} \frac{f(\tau)}{\sqrt{t_n - \tau}} d\tau\,.
$$
Apply the trapezoidal rule to the first term on the right-hand side.
Integrate the second term by parts on the small interval $[t_{n-1}, t_{n}]$,
then apply the trapezoidal rule. The derivative of the function at $t_{n-1}$
is generated as a result of the process. To evaluate its value,
there are multiple options, for instance, forward, backward and central
differences. We have tested the above-mentioned finite difference methods
on the cases where exact solutions are available; they appear to have the
same level of accuracy and computational speed; however, forward difference
is our choice as it is slightly more accurate. Combining the two parts,
 we obtain the  expression to approximate the integral
\begin{equation}
\label{eqn:FDhybrid}
I_n =\Delta t \sum_{i = 2}^{n-1} \frac{f_i}{\sqrt{t_n - t_i}}
+ \frac{\Delta t}{2} \Bigl(\frac{f_1}{\sqrt{t_n - t_1}}
- \frac{f_{n-1}}{\sqrt{t_n - t_{n-1}}}\Bigr)
+ \sqrt{\Delta t}\, (f_{n-1} + f_n) \,.
\end{equation}	
The scheme introduced here is a combination of two techniques, during which
forward finite differencing is applied; we will thus refer to it as the FD
Hybrid method in this paper. In Nevskii and Osiptsov's study of the effects of
unsteady and history forces in the gravity convection of suspensions
\cite{Nevskii2008}, a similar scheme has been mentioned.

\subsection{Predictor corrector method}
\label{subsec: m2}
Noticing that often the numerator of the integrand involves time-derivatives
of various terms, instead of what has been shown in \eqref{intHisVer01},
the integral is rewritten in the following way:
\begin{equation}
I(t) = \int_0^t  \frac{\dot{f}(\tau)}{\sqrt{t - \tau}} \,  d\tau \,.
\label{intHisVer02}
\end{equation}
On each subinterval $[t_{j -1}, t_j]$, $j =2, \dots,  n$, approximate
$\dot f(\tau)$ by forward difference at  $t_{j-1}$ (which may also be
regarded as a central difference approximation at the midpoint of the interval);
then the rest of the integral can be integrated analytically.
 When applied to the entire interval, a weighted sum is obtained to approximate
the integral, namely,
\begin{equation}\label{eqn:PC}
I_n =
\frac{2}{\Delta t} \sum_{j = 2}^{n} \bigl[ \sqrt{t_n - t_{j-1}}
- \sqrt{t_n - t_j}\,\bigr] \, (f_j - f_{j-1}) .
\end{equation}
In the process of evaluating the time-dependent function at the next step,
a pre\-dictor-corrector scheme based on this formula for evaluating the integral
can then be applied effectively to improve accuracy. In combination with the
time-stepping formula, this scheme will be referred to as the PC Method.
This turns out to be similar to the approach of discretizing the history
integral in \cite{reeks:1573}, further developed by Alexander
\cite{MR2069946}.

\section{Physical problems in one dimension}
\label{sec:PhyApp}

To test the above schemes, we apply both to various physical problems
for which exact analytical solutions are attainable so that we can compare
the numerical results against the exact solutions thereby validating the methods.
 By way of comparison, selected numerical methods from the existing literature
are also applied to the same problems.

\subsection{Sedimentation}
\label{subsec: Sedi}
The first problem we consider is sedimentation under gravity. For a sphere
initially released from rest in an infinite stagnant fluid, upon
nondimensionalizing the particle's velocity with the steady Stokes settling
velocity: $(2a^2 g \Delta \rho)/(9\mu)$, and scaling time with the viscous
diffusion time: $a^2/\nu$ (most physical parameters were introduced in
Section \ref{sec:Introduction}; additionally $\nu=\mu/\rho_f$ is the kinematic
viscosity and $\Delta \rho$ is the difference between the densities of the
particle and fluid), the dimensionless equation of motion is given by
\cite{clift2005bubbles}
\begin{equation}\label{11-15}
\begin{gathered}
 \frac{1}{B}\frac{dW}{dt}
=1 - W - \frac{1}{\sqrt\pi} \int_0^t \frac{\dot W(\tau)}{\sqrt{t - \tau}}d\tau,
\\
W(0) = 0,  \quad  \dfrac{dW}{dt}(0) = B\,,	
\end{gathered}
\end{equation}
where $W(t)$ is the ratio of the instantaneous velocity to the steady settling
velocity in Stokes flow, and $B$ is the dimensionless acceleration, $dW/dt$,
at $t = 0$, which depends solely on the density ratio $\gamma = \rho_p/\rho_f$,
specifically, $B = {9}/({2\gamma + 1})$. Note that the second initial
condition on $dW/dt$ can be inferred from the differential equation itself
by evaluating both sides at $t=0$ and using the first initial condition on $W$.
Thus, strictly speaking, it is not needed. Depending on whether the density
ratio $\gamma$ is less than, equal to, or greater than the critical value
$\gamma_c=\frac58$ (corresponding to $B_c=4$),  \eqref{11-15} has the
analytical solution
\[
W(t) =
\begin{cases}
1 - \dfrac{a}{a - b} \exp({b^2 t}) \,\text{Erfc}(b\sqrt{t}) - \dfrac{b}{b - a}
\exp({a^2 t})\, \text{Erfc}(a\sqrt{t})
& B>4, \\
1 + (8t-1)\exp(4t)\text{Erfc}(2\sqrt{t})-4\sqrt{t/\pi}
& B=4, \\
1 - \Re\{W[(Y+iX)\sqrt{t}]\} - \dfrac{X}{Y}\Im\{W[(Y+iX)\sqrt{t}]\}
& B<4,
\end{cases}
\]
where $a = \frac{B}{2} (1 + \sqrt{1 - {4}/{B}})$,
$b = \frac{B}{2}  (1 - \sqrt{1 - 4/B})$, $X=\frac{B}{2}$,
$Y=\frac{B}{2}\sqrt{4/B-1}$ and $W(Z)$ is the complex error function
$W(Z)=\exp(-Z^2) \bigl( 1+\frac{2 i}{\sqrt{\pi}} \int_0^Z \exp(\xi^2)\,d\xi\bigr)$,
with $\Re$ and $\Im$ referring to the real and imaginary parts.
 Erfc is the complementary error function defined by
$\text{Erfc}(z) \equiv 1 - \text{Erf}(z) = \frac{2}{\sqrt{\pi}} 
\int_z^{\infty} \exp(-t^2) dt $.
The expression for $B<4$ is also equivalent to the real part of the expression
given for $B>4$. Both reduce to the middle expression when $B$ approaches 4.


\subsubsection{Numerical solution}

To approximate $W(t)$ at the $(n+1)$-st time step, forward differencing is applied
to the derivative in \eqref{11-15} resulting in the update equation
\begin{equation}
w_{n+1} = B\Delta t + (1 - B\Delta t) w_n - \frac{B\Delta t}{\sqrt \pi} \, I_n\,,
\label{eq: Wn}
\end{equation}
where $w_n$ is the numerical solution at $t_n$ and $I_n$ is the numerical
approximation of the integral $\int_0^t \frac{\dot W(\tau)}{\sqrt{t - \tau}} d\tau$
at the $n$-th time step. We then apply  \eqref{eqn:FDhybrid}
and \eqref{eqn:PC} to the integral and obtain trajectories accordingly.

\subsubsection{Implementation of FD Hybrid and PC methods}
We now demonstrate how to implement the two methods in full detail for
 the 1-D sedimentation problem. Detailed descriptions for the other problems
will not be given, since they are quite similar to this one.
\smallskip

\noindent\textbf{FD Hybrid:} Apply \eqref{eqn:FDhybrid} to $I_n$ in \eqref{eq: Wn}.
Notice that \eqref{eqn:FDhybrid} requires $n \ge 3$, that is, if we denote
the discretized velocity as the sequence $w_i$ ($i = 1, 2, \dots$),
we need to know $w_1, w_2, w_3$ before we can apply the update equation.
 From the initial conditions, it is natural to let $w_1 = 0$,
 $w_2 = B\Delta t + w_1$ by first-order forward finite differencing, and
 $w_3 = -2\big(B \Delta t + \frac{3}{2}w_1 - 2w_2\big)$ by the second-order
forward finite differencing for the first derivative. Effectively, this
 ends up approximating the solution as a straight line along the first three
time nodes. Then $I_3$ is given by
\begin{align*}
I_3 & =
\Delta t  \frac{\dot w_2}{\sqrt{t_3 - t_2}}
+ \frac{\Delta t}{2} \Bigl(\frac{\dot w_1}{\sqrt{t_3 - t_1}}
- \frac{\dot w_{2}}{\sqrt{t_3 - t_{2}}}\Bigr)
+ \sqrt{\Delta t}\, (\dot w_{2} + \dot w_3)
\\
& = \frac{w_3 - w_2}{\sqrt{\Delta t}}
+ \frac{1}{2}\Bigl(B\sqrt{\frac{\Delta t}{2}}
- \frac{w_3 - w_2}{\sqrt{\Delta t}}\Bigr)
+ 2\frac{w_3 - w_2}{\sqrt{\Delta t}}
\\
& = \frac{B\sqrt{\Delta t}}{2\sqrt 2}+ \frac{3(w_3 - w_2)}{2\sqrt{\Delta t}}\,.
\end{align*}
Here the derivatives are approximated by forward difference, except that we use
 backward on the term $\dot w_3$.
In general, for $n \ge 3$, the updating expression for the memory integral is
\begin{align*}
I_n
& =
\Delta t \sum_{i = 2}^{n-1} \frac{\dot w_i}{\sqrt{t_n - t_i}}
+ \frac{\Delta t}{2} \Bigl(\frac{\dot w_1}{\sqrt{t_n - t_1}}
- \frac{\dot w_{n-1}}{\sqrt{t_n - t_{n-1}}}\Bigr)
+ \sqrt{\Delta t}\, (\dot w_{n-1} +\dot w_n)\\
& =
\sum_{i=2}^{n-1}\frac{w_{i+1} - w_i}{t_n-t_i}
+ \frac{1}{2} \Bigl(\frac{B\Delta t}{\sqrt{t_n -t_1}}
- \frac{w_n - w_{n-1}}{\sqrt{\Delta t}}\Bigr)
+  \frac{2(w_n - w_{n-1})}{\sqrt{\Delta t}}\,,
\end{align*}
which can be substituted into \eqref{eq: Wn} to complete the implementation.
\smallskip

\noindent\textbf{PC method:} Apply \eqref{eqn:PC} to $I_n$ in \eqref{eq: Wn}.
The PC Method requires only two initial points to get started as suggested in
\eqref{eqn:PC}, that is, $w_1 = 0$ and $w_2 = B\Delta t + w_1$.
Thus after substitution and simplification,
$$
I_2 =
\frac{2}{\sqrt{\Delta t}}  \, (w_2 - w_{1}) .
$$
In general for $n \ge 2$, the update expression to estimate the memory integral
in 1-D sedimentation is given by
\begin{equation}\label{eq: WnIn}
I_n =
\frac{2}{\Delta t} \sum_{j = 2}^{n} \big[ \sqrt{t_n - t_{j-1}}
- \sqrt{t_n - t_j}\big] \, (w_j - w_{j-1}).
\end{equation}
However, in the PC Method, two evaluations of the integral are required
for each time step. Suppose we have obtained the sequence $w_1, w_2, \dots, w_n$.
To approximate $w_{n+1}$, the first step is to obtain the approximation for the
history integral $I_n$ using \eqref{eq: WnIn}. In the ``predictor''
step, we continue by obtaining an approximation to $w_{n+1}$ by \eqref{eq: Wn};
let us call this $w_{n+1}^*$. We then use the resulting length $(n+1)$
sequence of $w_i$ to evaluate the integral again, denoted as $I_{n+1}^*$.
The ``predictor'' step then produces the final desired $w_{n+1}$ using
\[
w_{n+1}  = B\Delta t + (1 - B\Delta t) \frac{w_n+w_{n+1}^*}{2}
- \frac{B\Delta t}{\sqrt \pi} \, \frac{I_n + I_{n+1}^{*}}{2}.
\]

In addition to the above two methods, Brush's method and Daitche's first-order
scheme are also implemented for comparison. Figure \ref{fig:compSedi} shows the
 sedimentation velocity of a particle traveling until a dimensionless time of
 50 when $\gamma = 1/2 < 5/8$ with a uniform time step of 0.01.
The analytic solution is calculated and shown by the thick purple line.
The numerical results from different methods are differentiated by colors:
FD Hybrid (red), PC method (green), Brush's method (blue) and Daitche's first-order
scheme (yellow).

\begin{figure}[ht]
\begin{center}
\subfloat[$W(t), t \in (0, 50)$]{\label{fig:compSed}
\includegraphics[width=0.49\textwidth]{fig1a} % 20140503compSed01}
}
\subfloat[Error Plot]{\label{fig:errSed}
\includegraphics[width=0.49\textwidth]{fig1b} % 20140215_compSed_pic02}
}
\end{center}
\caption{1-D Sedimentation at $\gamma=0.5$ via FD Hybrid, PC,
Brush and Daitche's first-order methods.}
\label{fig:compSedi}
\end{figure}

As seen in Figure \ref{fig:compSed}, the numerical results all agree with
the exact solution. The particle is released from rest, then accelerates until
it reaches at a terminal velocity. The curves overlap with one another,
making it difficult to differentiate among them. The error plot in
Figure \ref{fig:errSed} helps highlight their differences. Here, the error
is defined as the difference between the exact and numerical solutions.
Among the four methods, PC is distinguished by having errors that are smaller
from the beginning and later become closest to zero, while Brush's method
strays farthest away from the exact solution at the very beginning.
 Both FD Hybrid and Daitche's first-order start off with smaller errors than
Brush's; however, for long enough times, all three seem to agree with each other,
 which suggests that the FD Hybrid method is of first-order accuracy.
We also observe that the FD Hybrid, Brush and Daitche's methods always
underestimate the solution, for their errors stay negative for all time.

In addition, we show the ``integrated error'', i.e., the integral of the magnitude
 of the instantaneous error over the entire time interval in
Table~\ref{tbl:IEsedi}.

\begin{table}[ht]
\caption{Integrated Errors (Sedimentation)} \label{tbl:IEsedi}
\begin{center}
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{|l|l|l|l|l|}
\hline
Integrated Error  & FD Hybrid & PC & Brush & Daitche
\\
\hline
$\Delta t =  0.01$ & 0.15713 & 0.00336 & 0.19070 & 0.16053
\\
\hline
$\Delta t =  0.001$ & 0.01555 & 0.00039 & 0.02453 & 0.01620
\\
\hline
\end{tabular}
\end{center}
\end{table}

We calculate these for all four numerical methods when $\Delta t = 0.01$ and
$0.001$.  It is not surprising that the integrated errors decrease as the
time step gets
finer for all cases. The PC method outperforms the other schemes.

\begin{figure}[ht]
\begin{center}
 \includegraphics[width=0.6\textwidth]{fig2} % 20140506compSedLog01
\end{center}
\caption{Log-Log plot of the integrated error versus $\Delta t$.}
\label{fig:loglog}
\end{figure}

The Log-Log plot of the errors versus $\Delta t$ shown in
 Figure~\ref{fig:loglog} displays the order of accuracy clearly.
The PC method is second-order accurate, while FD Hybrid, Daitche,
and Brush are first order.

\subsection{Impulsive motion}
\label{subsec: Imp}

The second test problem models the case where the fluid is assumed to be
 at rest and the particle is given an initial velocity and allowed to
relax freely in the viscous fluid, in the absence of gravity.
The motion of the spherical particle is in one direction (e.g., parallel
 to unit vector {\bf i} in Cartesian coordinates) and buoyancy effects are
neglected. As we shall see below, the initial naive formulation of this
 problem leads to unphysical results, although it is mathematically well-posed
and can still be used for benchmarking of the numerical results.
This unphysical result turns out to be due to the infinite acceleration that
would be implied by the mismatch between the fluid and particle velocities
at the initial instant. The resolution of this paradox, shown in the second
formulation, involves replacing the initial velocity of the particle with an
impulsive force that acts on it, causing it to attain the desired velocity
over a very short (nearly zero) time span.

\subsubsection{Formulation I}
Substitute $\mathbf{u} = 0$ (for a quiescent fluid) into  \eqref{eqn:eom};
the equation of motion then reads
\begin{equation}
\begin{gathered}
 m_p \dot{\mathbf{V}\,}
=-\frac{m_f}{2} \dot{\mathbf{V}\,}
- 6\pi \mu a \mathbf{V}
- 6a^2 \sqrt{\pi \mu \rho_f} \int_0^t \frac{\dot {\mathbf{V}\,}(\tau)}{\sqrt{t - \tau}} \, d\tau
\\
 {\mathbf{V}}(0) = v_o \,{\bf i},
\end{gathered}
\label{eqn:yEqn}
\end{equation}
where $\mathbf{V}$ is the velocity of the particle, and $v_o$ is its initial
velocity. Let $\mathbf{V}(t) = {\bf i}\, y(t)$,
choose the velocity scale $v_0$ and the time scale $(m_p + m_f /2)/(6 \pi \mu a)$
to render $y$ and $t$ dimensionless,
substitute into  \eqref{eqn:yEqn} and rename the dimensionless variables $y$
and $t$ again for simplicity, to get
\begin{equation}
\begin{gathered}
\dot y = -y - \beta \int_{0}^t \frac{\dot y(\tau)}{\sqrt{t - \tau}} d\tau
\\
y(0)  = 1.
\end{gathered} \label{eqn:yEqnDimLess}
\end{equation}
The dynamics are determined by a single dimensionless parameter
\[
\beta = \sqrt{\frac{9\rho_f}{\pi(2\rho_p + \rho_f) }}
\]
that solely depends on the particle and fluid densities and whose
value ranges from 0 to $3/\sqrt{\pi}$, as the ratio $\rho_f/\rho_p$
goes from zero to infinity.
\smallskip		

\noindent\textbf{Analytic solution:}
Laplace transforms lead to the analytic solution (see Section \ref{sec:appdFormI}
for details)
\begin{equation}\label{eqn:exSolYfor1}
y(t)=  \sum_{i = 1}^2 A_iX_i\exp (X_i^2 t) \text{Erfc}\bigl( - X_i\sqrt{t}\,\bigr).
\end{equation}
Erfc is the complementary error function defined earlier. Coefficients
$A_i$ and arguments $X_i$ are given by \eqref{X012} and \eqref{A012} in the appendix.
They are defined in terms of the parameter $b = \beta\sqrt{\pi}/2$.
To better understand this solution, let us examine its asymptotic behavior
for both short and long times.
\smallskip		

\noindent\textbf{Short-time behavior (near $0^+$):}
For each term in \eqref{eqn:exSolYfor1}, the series expansion for small $t$
looks like
\begin{equation}
\exp\bigl(X^2\, t\bigr) \text{Erfc} \bigl(-X \, \sqrt{t}\,\bigr)
= 1 + \frac{2X\sqrt{t}}{\sqrt{\pi}} + X^2 t + \frac{4X^3 t^{3/2}}{3\sqrt{\pi}}
+ \mathcal O\bigl(t^2\bigr).
\label{eqn:nearzero}
\end{equation}
Thus,
\begin{align*}
y(t)
& = (A_1X_1 + A_2 X_2)
+ \sqrt{t}\frac{2}{\sqrt{\pi}}\bigl(A_1X_1^2 + A_2 X_2^2\bigr)
\\
& \quad + t \bigl(A_1X_1^3 + A_2 X_2^3\bigr)
+ t^{3/2} \frac{4}{3\sqrt{\pi}} \bigl(A_1X_1^4 + A_2 X_2^4\bigr)
+ \mathcal O \bigl(t^2 \bigr)
\\
& = 1 - t + \mathcal O\bigl(t^{3/2}\bigr).
\end{align*}
Here, we have used the results
$\sum_{i = 1}^2 A_i X_i = 1$,
$\sum_{i = 1}^2 A_i X_i^2 = 0$,
$\sum_{i = 1}^2 A_i X_i^3 = -1$,
$\sum_{i = 1}^2 A_i X_i^4 = 2b$.
\smallskip		

\noindent\textbf{Long-time behavior (near $\infty$):}
The large-$t$ series expansion of each term has the form
\begin{equation}
\exp\bigl(X^2 t\bigr) \text{Erfc} \bigl(-X \sqrt{t}\bigr)
= \frac{1}{X} \Bigl( - \sqrt{\frac{1}{t\pi}} + \frac{1}{2\sqrt \pi X^2} t^{-3/2}
+ \mathcal O\bigl(t^{-5/2}\bigr)\!\Bigr).
\label{eqn:nearinf}
\end{equation}
As such,
\begin{align*}
y(t)
& =- t^{-1/2} \sqrt{\frac{1}{\pi}} (A_1 + A_2)
+ t^{-3/2} \frac{1}{2\sqrt\pi} \Bigl(\frac{A_1 }{X_1^2 } + \frac{A_2}{X_2^2}\Bigr)
+ \mathcal O\bigl(t^{-5/2}\bigr)
\\
& = 2b\sqrt{\frac{1}{\pi}} t^{-1/2} + \mathcal O\bigl(t^{-3/2}\bigr)
\\
& = \frac{\beta}{\sqrt{t}} + \mathcal O\bigl(t^{-3/2}\bigr)
\end{align*}
since $A_1 + A_2 = - 2b$.

The fact that the longtime asymptotic form of the particle velocity decays
relatively slowly as $t^{-1/2}$ leads to an unphysical result.
Namely, if we consider the particle displacement, which is the integral
of velocity with respect to time, the resulting integral diverges,
which suggests that the particle will travel an infinite distance! Therefore,
while the above impulsive motion problem is mathematically acceptable and
 useful for benchmarking the numerical methods, it does not correctly model
the physics of particle motion. To resolve this paradoxical issue, we
re-examine this problem below using a formulation that takes into 
account the force of the impulse delivered to the particle.
\smallskip		

\noindent\textbf{Relationship to the 1-D sedimentation:}
Before we introduce the second formulation of the impulsive motion example,
let us point out the strong similarity between this problem and the previous 1-D
sedimentation. Consider  \eqref{11-15} and \eqref{eqn:yEqnDimLess}
and notice their resemblance.
The two dimensionless groups are themselves related by: $\beta = \sqrt{{B}/{\pi}}$.
If we rescale time by $B$ in the sedimentation problem, i.e., $T = B  t$,
let $\tau'$ denote the scaled dummy variable of time-integration,
and further introduce a variable $w(t) = 1- W(t)$, \eqref{11-15} becomes
\begin{equation}
\frac{dw}{dT} =  -w - \sqrt{\frac{B}{\pi}} \int_0^{T}
\frac{\dot w(\tau')}{\sqrt{T - \tau'}} \,d\tau' .
\label{eqn:eqv}
\end{equation}
Equation \eqref{eqn:eqv} appears to be equivalent to \eqref{eqn:yEqnDimLess}.
To confirm that they are mathematically the same, we also need to examine the
initial conditions. Since $w(0)=1-w(0)=1$,
the initial condition on $w$ turns out to be the same as $y(0) = 1$.
As for the initial acceleration,
$$
\frac{dw}{dT}\bigr|_{t = 0}
= \frac{d(1 - W(t))}{d(Bt)}\bigr|_{t = 0}
= -\frac{1}{B}\frac{dW}{dt}\bigr|_{t = 0} = -1\,.
$$
The initial accelerations are thus also consistent because
\[
\dot y(0)
= - y(0) - \beta\int_0^0 \frac{\dot y(\tau)}{\sqrt{t - \tau}} d\tau
=  -1\,.
\]
We see that the two problems are actually equivalent to each other mathematically,
which raises an interesting question.
Given the unphysical result from the first impulsive motion formulation,
does a similar issue affect the sedimentation problem? The above equivalency
seems to lead to the conclusion that while in the sedimentation problem,
the particle eventually does relax to its final constant Stokes settling velocity,
if the distance that it would lag behind another particle that always settled
at that terminal velocity were measured, that distance would diverge in time.
In other words, in an infinite container, if one particle is settling at the
Stokes settling velocity, and another one is released from rest as the first
particle passes by, the second particle will lag behind and the distance between
them would diverge as $t\rightarrow\infty$. Therefore, there is a somewhat
 unphysical aspect to the sedimentation problem with the history force as well.
This subtle issue appears not to have been noticed previously.
\smallskip		

\noindent\textbf{Numerical results:}
In the exact solution for $y(t)$, it appears that $\beta$ has a critical value
at $2/\sqrt{\pi}$ (where the quantities under the square roots undergo changes
of sign). However, Coimbra and Rangel \cite{MR1647206} pointed out 
that the critical value is
only of mathematical relevance and does not imply a change of physical character
of the problem. For the numerical experiments, we set $\beta = 1/\sqrt{\pi}$
 which means that the particle is four times as dense as the fluid.
Once again, four numerical methods are implemented: FD Hybrid, PC,
Brush and Daitche's first-order methods. The scaled velocities computed
by each method are plotted in Figure~\ref{fig:compY}.

\begin{figure}[ht]
\begin{center}
\subfloat[$y(t), t \in (0, 50)$]{\label{fig:compy}
\includegraphics[width=.49\textwidth]{fig3a} % 20140215_compY_pic01
}
\subfloat[Error Plot]{\label{fig:errY}
\includegraphics[width=.49\textwidth]{fig3b} % 20140215_compY_pic02
}
\end{center}
\caption{Impulsive Motion (Formulation I): FD vs.\ PC vs.\ Brush vs.\
Daitche's first-order method.}
\label{fig:compY}
\end{figure}

The particle slows down substantially right after the impulsive start,
then approaches the expected time dependence of its velocity,
$\beta/\sqrt{t}$, asymptotically as time increases.
	
The PC method appears to yield the best result, as seen in
Figure \ref{fig:errY} which tracks the errors versus time.
The other three methods differ at the beginning; however, given long enough time,
their errors tend to approximately the same level. The integrated errors
are displayed in Table~\ref{tbl:IEy} for two choices of time step.	

\begin{table}[ht]
\caption{Integrated Errors (Impulsive Motion: Formulation I)}
\label{tbl:IEy}
\begin{center}
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{|l|l|l|l|l|}
\hline
Integrated Error  & FD Hybrid & PC & Brush & Daitche
\\
\hline
$\Delta t =  0.01$ & 0.03685 & 0.00418 & 0.06459 & 0.03863
\\
\hline
$\Delta t =  0.001$ & 0.00330 & 0.00042 & 0.01207 & 0.00389
\\
\hline
\end{tabular}
\end{center}
\end{table}
		
Besides accuracy, the execution time of a numerical method is another
important metric to evaluate its performance. For this purpose, we record
the total computational time for each method and use the lowest one (FD Hybrid)
as the unit of measurement.
We find that the FD Hybrid and Brush's method both take approximately one
time unit to run, while it takes the PC method about 1.5 times as long,
and Daitche's method requires almost 6 time units.

Both of the proposed methods appear to work as they produce results
that agree with the analytic solutions. On the merits of accuracy alone,
the PC Method has an obvious edge. As for computational cost, FD Hybrid seems
to be a little better. However, if the overall execution time is not prohibitively
long, one can sacrifice the computation cost to achieve the better accuracy. 	 	

\subsubsection{Formulation II}
In \cite{maxey:883}, it is assumed that the initial velocities of the particle
and the fluid are the same. But in the impulsive motion problem, while the
fluid is at rest, the particle is given an initial velocity $v_0 \neq 0$,
which suggests that a direct application of Maxey \& Riley's equation might
not be appropriate. As an alternative formulation, we can imagine that the
particle and fluid both start from rest, but that in addition to the Stokes
drag and the history force, an impulsive force of the form $S\,\delta (t-0^+)$
should be included in the right-hand side of the equation of motion.
The value $0^+$ in the argument of the delta function suggests that the
instantaneous impulse is applied ever so slightly after $t=0$.
This avoids potential confusion in applying the initial condition.
 In this form, the initial velocity of the particle coincides with that of the
fluid, so that we will be able to apply Maxey \& Riley's equation to get
\begin{equation}\label{eqn:yEqnNew}
\begin{gathered}
 (m_p+\frac{m_f}{2})\, \dot y = -{6\pi \mu a } \, y
- {6a^2 \sqrt{\pi \mu \rho_f}} \int_0^t \frac{\dot y(\tau)}{\sqrt{t - \tau}} d\tau
+ {S} \, \delta (t-0^+) ,
\\
 y(0)= 0.
 \end{gathered}
\end{equation}
Using the same time scale as before, ${(m_p + m_f/2)}/{6\pi \mu a}$, and the
new velocity scale $ S/(m_p+m_f/2)$,
we obtain the dimensionless equation of motion
\begin{equation}
\begin{gathered}
\dot y = -y - \beta \int_{0}^t \frac{\dot y(\tau)}{\sqrt{t - \tau}} d\tau
+ \delta (t-0^+),\\
y(0)  = 0.
\end{gathered} \label{eqn:yEqnNewDimLess}
\end{equation}
It should be noted that a delta function $\delta(t)$ has dimensions of
the reciprocal of its argument, in this case time. If time is scaled
like $t=T\hat{t}$, we have that $\delta(T\hat{t})=(1/|T|)\delta(\hat{t})$.
\smallskip

\noindent\textbf{Analytic solution:}
Using the Laplace transform, we can solve \eqref{eqn:yEqnNewDimLess} analytically:
\begin{equation}\label{eqn:exSolYFor2}
y(t) = \sum_{i = \pm}B_iX_i \exp(X_i^2 t)\text{Erfc} (-X_i \sqrt t)
\end{equation}
with the coefficients $B_i$ and factors $X_i$ given in \eqref{XpmBpm} in
Section \ref{sec:appdFormII}. A similar asymptotic analysis can be performed to
find the form of $y(t)$ when $t$ is close to zero and when it approaches
infinity, respectively.
\smallskip

\noindent\textbf{Short-time behavior (near $0^+$):}
Using \eqref{eqn:nearzero}, when $t \to 0^+$	
\begin{align*}
y(t)
& =\sum_{i = \pm} B_i X_i
+ t^{1/2}\frac{2}{\sqrt\pi} \sum_{i = \pm} B_i X_i^2
+ t \, \sum_{i = \pm} B_i X_i^3
+ t^{3/2} \, \frac{4}{3\sqrt \pi} \sum_{i = \pm} B_i X_i^4
+ \mathcal O(t^2) \\
& = 1 - 2\beta \sqrt{t} + \mathcal O(t).
\end{align*}
Here we have used
$\sum_{i = \pm} B_iX_i = 1$,
$\sum_{i = \pm} B_iX_i^2 = -2b$,
$\sum_{i = \pm} B_iX_i^3 = -1 + 4b^2$.
\smallskip

\noindent\textbf{Long-time behavior (near $\infty$):} By \eqref{eqn:nearinf},
when $t \to \infty$
\begin{align*}
y(t)
& =-(B_++B_-) \sqrt{\frac{1}{t\pi}}
+ \Bigl(\frac{B_+}{X_+^2} + \frac{B_-}{X_-^2}\Bigr)\cdot
\frac{1}{2\sqrt\pi}\cdot t^{-3/2} + \mathcal O\bigl(t^{-5/2}\bigr)
\\
& = \frac{\beta}{2} t^{-3/2}+ \mathcal O\bigl(t^{-5/2}\bigr).
\end{align*}
Note that $B_+ + B_- = 0$,  ${B_+}/{X_+^2} + {B_-}/{X_-^2} = \beta \sqrt\pi$.
\smallskip

\noindent\textbf{Displacement:}	
With the long-time form of the velocity now decaying in time like $t^{-3/2}$,
its integral with respect to time as the upper limit tends to infinity yields
a finite value. Thus the unphysical result associated with the previous
formulation does not plague this formulation of the impulsive start problem.
For the purpose of numerical validation, both formulations can be used since
they are mathematically well-posed problems. But the latter formulation is
the correct physical one for the impulsive start example.
\smallskip

\noindent\textbf{Numerical results:} Because of the presence of the Dirac delta
function $\delta (t-0^+)$, the above-mentioned numerical methods cannot
be applied directly until a scheme for dealing with the delta function
numerically has been devised. In order to make the system more amenable to
a numerical approach, first express the solution in the form
$y(t) = \alpha H(t-0^+)+ w(t)$ where $H(t)$ is the Heaviside step function,
defined by $ H(t) = \int_{-\infty}^ t \delta (s)\, ds$.
The derivative of the solution is thus
$\dot y (t) = \alpha \delta (t-0^+) + \dot w(t)$.
Substitute this into \eqref{eqn:yEqnNewDimLess} and balance the Delta
 functions on both sides of the equation by setting $\alpha =1$ to get,
for $t > 0^+$,
\begin{equation}\label{eqn:wt}
\dot w(t)  = -1- w(t) - \frac{\beta}{\sqrt t}
- \beta \int_{0}^t \frac{\dot w (\tau)}{\sqrt{t - \tau}} \,d\tau\,.
\end{equation}
The initial condition $y(0) = 0$ (applied at $t=0<0^+$) implies that $w(0) =0$.
This provides an explanation for the appearance of the terms
$(v_i-u_i)/\sqrt t$ as part of the history force in \cite{michaelides:1579}
and $w_i(0^+)/\sqrt{t}$ in \cite{Vojir1994547}. However, there is a new term
in the equation with singularity at $t = 0$. We choose the predictor-corrector
(PC) method to treat this problem because we can modify it to evaluate the
function at the midpoint of a time step so that it avoids the initial
singularity. Specifically, when the scheme is applied to
$\dot y = f(y, t) + r(t)$ where $r(t)$ has an explicit form in $t$
(in our case, $r(t) = -\beta/\sqrt t$ ), we modify the PC method to
\begin{align*}
y^\ast
& = y_n + \Delta t \bigl[f^n + r^{n+\frac{1}{2}}\bigr]
\\
y^{n+1}
& = y_n + \Delta t \bigl[\frac{f^n + f^{\ast}}{2} +  r^{n+\frac{1}{2}}\bigr]
\end{align*}
where $f^n = f(y_n, t_n)$, $f^\ast = f(y^*, t_n)$.

\begin{figure}[ht]
\begin{center}
\subfloat[Exact vs.\ Numerical (midPC)]{\label{fig:Form2traj}
\includegraphics[width=.5\textwidth]{fig4a} % 20140215_w_alone_pic01
}
\subfloat[Errors]{\label{fig:Form2err}
\includegraphics[width=.5\textwidth]{fig4b} % 20140215_w_alone_pic02
}
\end{center}
\caption{Exact and numerical solutions $y(t)$ and the corresponding error.}
\label{fig:impMo}
\end{figure}
	
The numerical results obtained with a time step of 0.01 are shown in
Figure \ref{fig:impMo}. The plot is obtained by solving  \eqref{eqn:wt}
numerically to obtain $w(t)$ and calculating the resulting $y(t)$
from the relation $y(t)=1+w(t)$ for $t>0$. The error is largest at the very
beginning as illustrated in Figure \ref{fig:Form2err}, which makes sense
since the abrupt change (due to the impulsive force) occurs when $t = 0$.
The integrated error in this case is 0.0273, and as the time step is
reduced to 0.001, the integrated error decreases to 0.0097, roughly
decreasing by a factor of $\sqrt{10}$.

\subsection{Oscillatory Sliding Motion}
\label{subsec: Osc}
Here we introduce another physical problem for which a more interesting
analytical solution can be obtained. Imagine that a spherical particle of
mass $m_p$ is placed in a cartridge filled with fluid which moves in a
horizontal sliding oscillation, in the absence of gravity (i.e., without
any sedimentation). The cartridge displacement is described by
$\Delta \sin (\omega t)\mathbf{i}$, where $\omega$ is the temporal frequency of
the linear oscillations and $\Delta$ denotes the amplitude.
We expect that the particle will eventually also undergo a back and
forth motion like its container, though not necessarily of the same
amplitude or in phase with the cartridge. We take the particle to have
the same initial velocity and acceleration as the fluid. Applying these to
\eqref{eqn:eom} and nondimensionalizing the system using time scale $1/\omega$
and length scale $\Delta$ result in the following dimensionless equation of motion:
\begin{equation}\label{eqn:eomX}
\begin{gathered}
 \ddot{X}(t) = - \beta \dot{X}(t) + ( 1- \alpha) \sin t
- \gamma \int_{0}^{ t} \frac{\ddot{X}(\tau) }{\sqrt{t - \tau}} \, d\tau,\\
X(0)= 0, \quad \dot X(0) = 0.
\end{gathered}
\end{equation}
$X(t)$ is the particle's displacement relative to the fluid
(i.e., in a frame of reference moving with the cartridge) and the
three dimensionless groups appearing in this equation are are defined as follows:
\[
\alpha =  \frac{3m_f}{2m_p + m_f}, \quad
\beta =  \frac{6\pi\mu a}{\omega(m_p + m_f/2)}, \quad
\gamma = \frac{6a^2 \sqrt{\pi \mu \rho_f/\omega} }{m_p + m_f/2}.
\]
Parameter $\alpha$ depends solely on the two densities, that is,
 $\alpha = 3\rho_f/(2\rho_p + \rho_f)$, and its value ranges from 0 to 3,
from extremely dense particles relative to the fluid to particles (or bubbles)
with very small densities; $\alpha=1$ corresponds to neutrally buoyant particles.
The second dimensionless group is the reciprocal of the Stokes number,
i.e., $\beta = {\text{St}^{-1}}$, which is a common parameter to characterize
the behavior of particles suspended in a fluid flow. It represents the ratio
of the relaxation time of the particle to the characteristic time of the flow:
the bigger the Stokes number, the longer it takes a particle to relax to
the same velocity as the fluid. The third parameter $\gamma$ quantifies
the importance of the history force. (Note that $\beta$ in this section is
 different from the one defined earlier.)
\smallskip

\noindent\textbf{Analytic solution:}
Once again Laplace transforms give the analytical solution (details are given
in Section \ref{sec:appdOSM}); namely,
\begin{equation}\label{eqn:exSolX}
X(t)  =  (1 - \alpha)\Bigl( \sum_{i = 1}^{6} A_i R_i \exp(R_i^2 t)
\text{Erfc}(-R_i\sqrt{t} ) + A_7  \Bigr),
\end{equation}
where
\begin{gather*}
R_1 = \exp\Bigl(\frac{\pi}{4} i\Bigr),\quad
R_2  = \exp\Bigl(\frac{3\pi}{4} i\Bigr),\quad
R_3 = \exp\Bigl(\frac{5\pi}{4} i\Bigr),\\
R_4 = \exp\Bigl(\frac{7\pi}{4} i\Bigr), \quad
R_{5,6} = \frac{-\gamma\sqrt\pi \pm \sqrt{\gamma^2 \pi - 4\beta}}{2}, \quad
R_{7,8} = 0.
\end{gather*}
along with the coefficients
$$
A_i = \frac{1}{ R_i^2\prod_{j=1,j\neq i}^{j=6} (R_i - R_j)} \quad (i=1,2,\dots, 6),
\quad
A_7 = \frac{1}{\beta} \,,
\quad A_8  =  \sum_{i=1}^ 6 \frac{1}{R_i}  \Bigm/{\prod_{i = 1}^{6} R_i} .
$$
\smallskip

\noindent\textbf{Numerical results:}	
FD Hybrid, PC, Brush's and Daitche's first-order methods are tested and compared
 with the analytic solution. The parameters are set to $\alpha = 0.3571$,
$\beta = 0.31$ and $\gamma = 0.3252$. The trajectories of the particle are
traced until a time of 50, or roughly 8 periods of oscillation of the container.

\begin{figure}[ht]
\begin{center}
\subfloat[ $X(t), t \in (0, 50)$]{\label{fig:}
\includegraphics[width=.5\textwidth]{fig5a} % 20140215_compX_pic01
}
\subfloat[Error Plot]{\label{fig:osErr}
\includegraphics[width=.5\textwidth]{fig5b} % 20140215_compX_pic02
}
\end{center}
\caption{Oscillatory Sliding Motion: FD vs.\ PC vs.\ Brush vs.\ Daitche's first-order}
\label{fig:oldSM}
\end{figure}
As shown in Figure~\ref{fig:oldSM}, all results agree and suggest that the
particle is initially ``thrown'' forward from its original position
by the starting motion of the container, but that after a while it stabilizes
and oscillates about a new equilibrium position. It is not surprising that
the errors oscillate as well as shown in Figure \ref{fig:osErr}.
From the error plot, Brush's method presents the largest errors,
while FD Hybrid and Daitche's method are of comparable magnitudes.
The PC method produces the smallest errors. This can be seen more clearly
in Table~\ref{tbl:IEX} which records the integrated errors when the
time step is 0.01 and 0.001 respectively.

\begin{table} [ht]
\caption{Integrated Errors (Oscillatory Sliding Motion)}
\label{tbl:IEX}
\begin{center}
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{|l|l|l|l|l|}
\hline
Integrated Error  & FD Hybrid & PC & Brush & Daitche
\\
\hline
$\Delta t =  0.01$ & 0.01740 & 0.00076 & 0.16238 & 0.01132
\\
\hline
$\Delta t =  0.001$ & 0.00375 & 2.33314e-05 & 0.05156 & 0.00015
\\
\hline
\end{tabular}
\end{center}
\end{table}


\section{Conclusions}
\label{sec:con}

Two methods were introduced to approximate history integrals of the form
$ \int_0^t {f(\tau)}/{\sqrt{t - \tau}} \, d\tau$ and
$ \int_0^t {\dot f(\tau)}/{\sqrt{t - \tau}} \, d\tau$ numerically.
 For the former integral, trapezoidal rule together with a forward difference
approximation of the derivative of $f$ (after integration-by-parts)
led to the FD Hybrid method; for the latter, the derivative $\dot{f}$
was treated by central differences at the mid-points of the time intervals,
leading to an approximation that could be used effectively in conjunction
with the Predictor-Corrector (PC) method to integrate the equations of
motion for particles affected by the history force. Both were applied
to three one-dimensional test problems that describe the physical problems
of sedimentation, impulsive motion and oscillatory sliding motion of a
particle in a fluid. The numerical results were compared with analytic
solutions obtained by Laplace transforms and both methods were shown to work
well in all cases.  We quantified their performance in terms of accuracy and
computation time, while comparing them to Brush's and Daitche's first-order methods
from the existing literature.

In all three cases, the PC method performed best, showing advantages in the
following respects: (i) the update expression is straightforward, making it easy
to implement; (ii) on the three one-dimensional problems tested, PC is the most
accurate; (iii) although it is not the fastest method, compared to the faster
FD Hybrid and Brush's methods, its computational time is within an acceptable
range; (iv) in the second formulation of the impulsive motion problem, it has
the advantage that it avoids having to evaluate the forcing function at its point
of singularity. Having validated these numerical schemes, they can be applied
with some confidence to more complex problems for which analytic solutions
are not available \cite{Xu2014}.

The test problems that we considered are useful for benchmarking the numerical
 methods in that their exact solutions could be obtained by Laplace transforms.
Our solution for the oscillatory sliding motion is new. For the impulsive start
problem, we provided two formulations, pointing out that the solution that
results by assuming that the particle has a different initial velocity than
the fluid, while mathematically acceptable, is not physically correct.
This pointed out a subtle issue with the interpretation of the existing
solution to the sedimentation problem as well, since it maps directly onto
the first formulation of the impulsive start example. The resolution of
the unphysical behavior involved introducing an impulsive force as a forcing
in the equation, while allowing the initial velocities of the particle and
the fluid to be the same.


\section{Appendix: Analytical solutions}\label{sec:appendix1}

We now provide the details for obtaining the exact solution for both
formulations of the impulsive motion example and the oscillating sliding
motion problem analytically by use of Laplace transform.

\subsection{Impulsive motion: formulation I}
\label{sec:appdFormI}

Apply Laplace transform to  \eqref{eqn:yEqnDimLess} and let $s$ denote
the Laplace transform variable. The Laplace transform of $y(t)$ is denoted
by $Y(s)$ and is found to be
$$
Y(s) =   \frac{ \sqrt{s} + \beta \sqrt{{\pi}}}{ s \sqrt{s}
+ \beta \sqrt{\pi }s +  \sqrt{s} }.
$$
The denominator is a cubic polynomial in $\sqrt{s}$ with the three roots:
\begin{equation}
\label{X012}
X_0 = 0, \quad
X_1 = {-b + \sqrt{b^2 - 1}}, \quad
X_2 = {-b - \sqrt{b^2 - 1}}\,,
\end{equation}
where $b = \beta\sqrt\pi/2$. $Y(s)$ can then be written as a partial fraction
expansion of the form
$$
Y(s) = \sum_{i = 0}^{2} \frac{A_i}{\sqrt{s} - X_i}
$$
in which the coefficients are given by
\begin{equation}
\label{A012}
A_0  = 2b, \quad
A_1  = \frac{1 - 2b^2 - 2b\sqrt{b^2 - 1}}{2\sqrt{b^2 - 1}}, \quad
A_2  = \frac{2b^2 - 1 - 2b \sqrt{b^2 - 1} }{2\sqrt{b^2 - 1}}.
\end{equation}
The inverse Laplace transform of $1/(\sqrt s - X_i)$ is 
$1/\sqrt{\pi t} + X_i\exp(X_i^2 t)\text{Erfc}\left(-X_i \sqrt{t}\right)$.
Noting that $ A_0 + A_1 + A_2 = 0$, the analytical solution for the
 impulsive motion problem simplifies to
\begin{align}\label{eqn:exSoly}
y(t) =    \sum_{i = 1}^2 A_iX_i\exp(X_i^2 t) \text{Erfc}\bigl( - X_i\sqrt{t}\,\bigr).
\end{align}


\subsection{Impulsive motion: formulation II}
\label{sec:appdFormII}

Upon taking the Laplace transform of  \eqref{eqn:yEqnNewDimLess}, we obtain
$$
sY = -Y - \beta s Y\sqrt{\frac{\pi}{s}} + 1 \quad \Rightarrow \quad
Y = \frac{1}{s + \beta \sqrt{\pi}\sqrt{s} + 1}.
$$
There are two roots (in $\sqrt{s}$) for the denominator, and to differentiate
them from the previous formulation, we denote them by $X_+$ \& $X_-$
with coefficients $B_+$ \& $B_-$ even though $X_{+(-)} \equiv X_{1(2)}$:
\begin{equation}
\label{XpmBpm}
X_{+,-} = -b \pm \sqrt{b^2 - 1}, \quad
B_{+,-} = \pm \frac{1}{2\sqrt{b^2 - 1}}.
\end{equation}
The partial fraction expansion of $Y(s)$ then reads
\[
 Y(s) = \frac{B_+}{\sqrt s - X_+} + \frac{B_-}{\sqrt s - X_-},
\]
whose inverse Laplace transform provides the exact solution given in
\eqref{eqn:yEqnNewDimLess}.

\subsection{Oscillatory sliding motion}
\label{sec:appdOSM}

Applying Laplace transform to \eqref{eqn:eomX},
\[
s^2{X(s)} = -\beta s{X} + \frac{1 - \alpha}{s^2 + 1}
- \gamma s^2{X}\sqrt{\frac{\pi}{s}}\,,
\]
where ${X} (s) = \mathcal L\{X(t)\}$. Solving for ${X}(s)$ we get
\[
{X}(s) = \frac{1- \alpha}{(s^2 + 1)(s^2 + \gamma\sqrt{\pi}s\sqrt{ s}+ \beta s ) }\,.
\]
Let $R = \sqrt s$ and find the roots of
$ {(R^4 + 1)(R^2 + \gamma\sqrt\pi R + \beta  )R^2}$. Denote these by $R_i$
and the corresponding coefficients in the partial fraction expansion of
$X(s)$ by $A_i$ $(i = 1, \dots, 8)$.
They are given explicitly following \eqref{eqn:exSolX}.
In terms of these, ${X}(s)$ can be written as
$$	{X}(s)
	= (1- \alpha)\Bigl( \sum_{i = 1}^{6} \frac{A_i}{\sqrt {s} - R_i}
+ A_7\frac{1}{{s}} + A_8 \frac{1}{\sqrt{s}}\Bigr).
$$
Noting that the inverse Laplace transform of ${1}/({\sqrt s - R_i })$ is
\[
 ( {1}/{\sqrt{\pi t}} + R_i \exp(R_i^2 t)\text{Erfc}(-R_i\sqrt{t} )),
\]
 that of
${1}/{\sqrt{s}}$ is ${1}/{\sqrt{\pi t}}$, and that of $ {1}/{s}$ is the
Heaviside step function, the solution for $X(t)$ is given by
$$	
X(t)  = (1 - \alpha)\Bigl( \sum_{i = 1}^{6} A_i\Bigl( \frac{1}{\sqrt{\pi t}}
+ R_i \exp(R_i^2 t)\text{Erfc}(-R_i\sqrt{t})\!\Bigr)
	+ A_7 H(t) + A_8 \frac{1}{\sqrt{\pi t}}\Bigr).
$$
Using the fact that $A_8 + \sum_{i = 1}^6 A_i  = 0$, this simplifies for
 positive time to the result given in \eqref{eqn:exSolX}.

\subsection*{Acknowledgement}
Shujing Xu thanks the Institute of Mathematical Sciences at Claremont
Graduate University for the Director's Fellowship that partially supported
her graduate studies.


\begin{thebibliography}{00}

\bibitem{MR2069946}
P.~Alexander;
\newblock {High order computation of the history term in the equation
of motion
  for a spherical particle in a fluid}.
\newblock { J. Sci. Comput.}, 21(2):129--143, 2004.

\bibitem{MR0134559}
A.~B. Basset;
\newblock {A treatise on hydrodynamics with numerous examples}.
\newblock Dover Publications Inc., New York, 1961.

\bibitem{Belmonte2000}
A.\ Belmonte, J.\ Jacobsen, A~Jayaraman;
\newblock Monotone solutions of a nonautonomous differential equation for a
  sedimenting sphere.
\newblock {Electronic Journal of Differential Equations}, 2001(62):1--17,
  2001.

\bibitem{Bombardelli2008}
F.\ Bombardelli, A.\ Gonz\'alez, Y.~Ni\~no;
\newblock Computation of the particle {B}asset force with a
  fractional-derivative approach.
\newblock {Journal of Hydraulic Engineering}, 134(10):1513--1520, 2008.

\bibitem{Boussinesq1885}
J.~V. Boussinesq;
\newblock Sur la r\'esistance qu'oppose un fluide ind\'efini au repos, sans
  pesanteur, au mouvement vari\'e d'une sph\`ere solide qu'il mouille sur toute
  sa surface, quand les vitesses restent bien continues et assez faibles pour
  que leurs carr\'es et produits soient n\'egligeables.
\newblock {Comptes Rendu de l'Academie des Sciences}, 100, 1885.

\bibitem{Brush1964}
L.~M.\ Brush, H.~W.\ Wo, B.~C. Yen;
\newblock Accelerated motion of a sphere in a viscous fluid.
\newblock {Journal of the Hydraulics Division}, 90:149--160, 1964.

\bibitem{candelier:1765}
F.\ Candelier, J.~R.\ Angilella, M.~Souhar;
\newblock On the effect of the {B}oussinesq--{B}asset force on the radial
  migration of a {S}tokes particle in a vortex.
\newblock {Physics of Fluids}, 16(5):1765--1776, 2004.

\bibitem{FLM:353169}
E.~J. Chang, M.~R. Maxey;
\newblock Unsteady flow about a sphere at low to moderate {R}eynolds number.
  part 1. oscillatory motion.
\newblock { Journal of Fluid Mechanics}, 277:347--379, 9 1994.

\bibitem{FLM:340128}
E.~J. Chang, M.~R. Maxey;
\newblock Unsteady flow about a sphere at low to moderate {R}eynolds number.
  part 2. accelerated motion.
\newblock { Journal of Fluid Mechanics}, 303:133--153, 10 1995.

\bibitem{clift2005bubbles}
R.\ Clift, J.~R.\ Grace, M.~E. Weber;
\newblock { Bubbles, drops, and particles}.
\newblock Dover Publications, 2005.

\bibitem{MR1647206}
C.~F.~M.\ Coimbra and R.~H. Rangel.
\newblock General solution of the particle momentum equation in unsteady
  {S}tokes flows.
\newblock {J. Fluid Mech.}, 370:53--72, 1998.

\bibitem{Daitche2012}
A.~Daitche;
\newblock Advection of inertial particles in the presence of the history force:
  Higher order numerical schemes.
\newblock Technical Report arXiv:1210.2576, 2012.

\bibitem{Daitche2011}
A.\ Daitche and T.~T\'el.
\newblock Memory effects are relevant for chaotic advection of inertial
  particles.
\newblock { Phys. Rev. Lett.}, 107:244501, 2011.

\bibitem{FLM:14341}
I.\ Kim, S.\ Elghobashi, W.~A. Sirignano;
\newblock On the equation for spherical-particle motion: effect of {R}eynolds
  and acceleration numbers.
\newblock { Journal of Fluid Mechanics}, 367:221--253, 6 1998.

\bibitem{ling:113301}
Y.\ Ling, J.~L.\ Wagner, S.~J.\ Beresh, S.~P.\ Kearney, and S.~Balachandar;
\newblock Interaction of a planar shock wave with a dense particle curtain:
  Modeling and experiments.
\newblock {Physics of Fluids}, 24(11):113301, 2012

\bibitem{LovBrady1993}
P.~M.\ Lovalenti and J.~F. Brady;
\newblock The hydrodynamic force on a rigid particle undergoing arbitrary
  time-dependent motion at small {R}eynolds number.
\newblock {Journal of Fluid Mechanics}, 256:561--605, 1993.

\bibitem{lukerchenko2010basset}
N.~Lukerchenko;
\newblock {B}asset history force for the bed load sediment transport.
\newblock In {Proceedings of the First IAHR European Division Congress,
  4--6 May 2010}, 2010.

\bibitem{maxey:883}
M.~R.\ Maxey and J.~J.\ Riley;
\newblock Equation of motion for a small rigid sphere in a nonuniform flow.
\newblock {Physics of Fluids}, 26(4):883--889, 1983.

\bibitem{CambridgeJournals:397249}
R.\ Mei, R.~J.\ Adrian, T.~J. Hanratty;
\newblock Particle dispersion in isotropic turbulence under {S}tokes drag and
  {B}asset force with gravitational settling.
\newblock {Journal of Fluid Mechanics}, 225:481--495, 1991.

\bibitem{michaelides:1579}
E.~E. Michaelides;
\newblock A novel way of computing the {B}asset term in unsteady multiphase
  flow computations.
\newblock {Physics of Fluids A: Fluid Dynamics}, 4(7):1579--1582, 1992.

\bibitem{citeulike:11142061}
E.~E. Michaelides;
\newblock Review - the transient equation of motion for particles, bubbles, and
  droplets.
\newblock {Journal of Fluids Engineering}, 119(2):233+, 1997.

\bibitem{Nevskii2008}
Y.~A.\ Nevskii, A.~N. Osiptsov;
\newblock The effect of unsteady and history forces in the gravity convection
  of suspensions.
\newblock {Moscow University Mechanics Bulletin}, 63:85--88, 2008.

\bibitem{reeks:1573}
M.~W.\ Reeks, S.~McKee;
\newblock The dispersive effects of {B}asset history forces on particle motion
  in a turbulent flow.
\newblock {Physics of Fluids}, 27(7):1573--1582, 1984.

\bibitem{hinsberg}
M.~A.~T.\ van Hinsberg, J.~H.~M.\ ten Thije~Boonkkamp, H.~J.~H. Clercx;
\newblock An efficient, second order method for the approximation of the
  {B}asset history force.
\newblock {J. Comput. Physics},  230:1465-1478, 2011.

\bibitem{visitskii2006}
E.~V.\ Visitskii, A.~G. Petrov;
\newblock Structurization of solid particles in a liquid medium under the
  action of a standing ultrasonic wave field.
\newblock {Doklady Physics}, 51:328--333, 2006.

\bibitem{Visitskii2009548}
Y.~V.\ Visitskii, A.~G.\ Petrov, M.~M. Shunderyuk;
\newblock The motion of a particle in a viscous fluid under gravity, vibration
  and {B}asset's force.
\newblock {Journal of Applied Mathematics and Mechanics}, 73(5):548 -- 557,
  2009.

\bibitem{Vodopaoyanov2010}
I.~S.\ Vodop'yanov, A.~G.\ Petrov, M.~M. Shunderyuk;
\newblock Unsteady sedimentation of a spherical solid particle in a viscous
  fluid.
\newblock {Fluid Dynamics}, 45:254--263, 2010.

\bibitem{Vojir1994547}
D.~J.\ Vojir, E.~E. Michaelides;
\newblock Effect of the history term on the motion of rigid spheres in a
  viscous fluid.
\newblock {International Journal of Multiphase Flow}, 20(3):547 -- 556,
  1994.

\bibitem{Xu2014}
S.~Xu;
\newblock {Effects of History and Lift Forces on Particle Trajectories in
  Oscillatory Rotating Fluids}.
\newblock PhD thesis, Claremont Graduate University, August 2014.

\bibitem{yang:1822}
S.~M.\ Yang, L.~G. Leal;
\newblock A note on memory-integral contributions to the force on an
  accelerating spherical drop at low {R}eynolds number.
\newblock {Physics of Fluids A: Fluid Dynamics}, 3(7):1822--1824, 1991.

\bibitem{yih1977fluid}
C.~S. Yih;
\newblock {Fluid Mechanics: A Concise Introduction to the Theory}.
\newblock West River Press, 1977.

\end{thebibliography}



\end{document}