Use of the Basset coefficient in the calculation of the velocity of a spheroid slowing down in a viscous incompressible fluid after a sudden impulse

The motion of a spheroid in a viscous incompressible fluid, after a sudden small impulse directed along the symmetry axis, may in principle be calculated from the linearized Navier-Stokes equations for fluid velocity and pressure, in combination with a boundary condition for the flow velocity at the surface of the spheroid. In the following we assume that the no-slip boundary condition holds. For a small body the no-slip condition must be relaxed and replaced by a partial slip boundary condition 1 .

The added mass tensor, which together with the body mass determines the initial velocity after the impulse, has long been known from the general expressions for an ellipsoid 2 . Similarly, the steady state friction tensor, which determines the steady state velocity for a constant force, is also known for the general ellipsoid 2 . The complete time-dependence of the velocity after a sudden impulse follows in principle by Fourier transform from the frequency-dependence of the friction tensor, but this is not known for a general ellipsoid. For axial motion of a spheroid the frequency-dependent friction coefficient can be found in principle from a set of coupled equations for the multipole moments of the force density induced on the spheroid, as shown first by Lai and Mockros 3 . It was shown by Lawrence and Weinbaum 4 , that the equations can be solved efficiently by truncation and matrix inversion.

At low frequency the friction coefficient for axial motion is linear in the square root of frequency, with a coefficient which can be found from the friction coefficient at zero frequency 5 ^,6 . As a consequence of this relation the velocity relaxation function shows a $t^{-3/2}$ long-time tail, independent of shape, size, and mass density of the spheroid. The transient effects on the short and intermediate time scale are less well known.

For the case of a sphere, Stokes 7 derived a simple explicit formula for the frequency-dependent friction coefficient. The impedance of a sphere is quadratic in the square root of frequency. The term linear in the square root is relevant both at low and at high frequency. The corresponding motion of a sphere has been studied by Boussinesq 8 . For a spheroid, the asymptotic behavior of the impedance at high frequency is given by a term linear in frequency, with coefficient given by the sum of spheroid mass and added fluid mass, and a term proportional to the square root of frequency, with a Basset coefficient. Lawrence and Weinbaum 4 calculated the latter from an expression derived by Landau and Lifshitz 9 and Batchelor 10 for a body of arbitrary shape. The Basset coefficient differs from the coefficient of the term linear in the square root of frequency at low frequency. In previous work we have calculated the Basset coefficient for a general ellipsoid for the three principal directions 11 .

The steady state friction coefficient and the effective mass are sufficient for the construction of a simple, approximate expression for the complex admittance, showing two poles in the complex plane of the square root of frequency. On the basis of their numerical work Lawrence and Weinbaum suggested a three-pole expression for the complex admittance corresponding to the axial motion of a spheroid, involving the steady-state friction coefficient, the effective mass, and the Basset coefficient 4 . It has been suggested that the approximate expression provides an adequate approximation to the actual admittance 12 ^,13 . In the following we examine this idea by comparison with a calculation of the admittance on the basis of the method of Padé approximants 14 . Our conclusion is that the simple two-pole expression is quite useful, and that the small deviations from the two-pole approximation are not well described by the Lawrence-Weinbaum expression.

We consider a prolate or oblate spheroid of mass $m_{p}$ immersed in a viscous incompressible fluid of shear viscosity $\eta$, mass density $\rho$. The fluid is assumed to fill the whole outer space, and a no-slip boundary condition holds at the surface of the spheroid. The spheroid and the fluid are at rest for $t<0$. We choose a cylindrical frame of coordinates $(r,\varphi,z)$ such that the center of the spheroid at rest is at the origin, and the $z$ axis is the axis of symmetry. Furthermore $r$ is the radial distance from the axis, and $\varphi$ is the azimuthal angle. We follow the convention of Happel amd Brenner 15 , and denote the longest semi-axis by $a$ and the shortest by $b$ in both the prolate and the oblate case. At time $t=0$ a sudden impulse $\mbox{\boldmath$S$}=S\mbox{\boldmath$e$}_{z}$ in the direction of the symmetry axis is imparted to the spheroid, corresponding to the applied force

The impulse is assumed to be sufficiently small that the resulting flow velocity $\mbox{\boldmath$v$}(\mbox{\boldmath$r$},t)$ and pressure $p(\mbox{\boldmath$r$},t)$ can be assumed to satisfy the linearized Navier-Stokes equations

The flow velocity $v$ tends to zero at infinity, and the pressure $p$ tends to a constant $p_{st}$. We are interested in the corresponding velocity $\mbox{\boldmath$U$}(t)=U(t)\mbox{\boldmath$e$}_{z}$ of the spheroid.

The equation of motion for the spheroid reads

where $K(t)$ is the hydrodynamic force exerted by the fluid, which may be calculated from the fluid stress tensor. In Fourier transform the equation becomes

Since $K_{\omega}$ is linear in $U_{\omega}$ this may be resolved to the form

with admittance $\mathcal{Y}(\omega)$. The latter may be cast in the form

Here $m_{a}$ is the added mass, which follows from the high-frequency behavior. The friction coefficient $\zeta(\omega)$ is at most linear in $\sqrt{\omega}$ at high frequency. The impedance is the inverse of the admittance.

The added mass $m_{a}$ follows from the theory of potential flow, and is well known 2 . The friction coefficient at zero frequency follows from the solution of the steady state Stokes equations. It takes the form

where $R$ has the dimension length. Its explicit expression is also well known 2 .

For a sphere of radius $a$ the friction coefficient takes the simple form 7

where $\alpha=\sqrt{-i\omega\rho/\eta},\;\mathrm{Re}\;\alpha>0$. Similarly, for a spheroid the friction coefficient can be expanded in powers of $\alpha$ as

The term linear in $\alpha$ can be expressed in terms of the coefficient at zero frequency $\zeta(0)=\zeta^{(0)}$ as

A similar theorem in tensor form holds for a body of arbitrary shape 5 ^,6 .

The first term in the asymptotic expansion of $\zeta(\omega)$ in powers of $\alpha^{-1}$ is linear in $\alpha$,

where $B$ is the dimensionless Basset coefficient 16 ^,4 . For a sphere $B=1$ from Eq. (8). The explicit expression for an ellipsoid can be evaluated from the theory of potential flow 2 by use of an argument due to Landau and Lifshitz 9 and to Batchelor 10 . The expression for an oblate spheroid was derived by Lawrence and Weinbaum 4 . It agrees with the one derived for a general ellipsoid 11 . The corresponding expression for a prolate spheroid has also been derived 11 ^,12 .

It is convenient to write the admittance in the form

where $\mathcal{Y}_{0}(\omega)$ is related to the Fourier transform of the hydrodynamic force by

It clearly is independent of the mass $m_{p}$.

We define the viscous relaxation time $\tau_{v}=a^{2}\rho/\eta$, and the corresponding complex variable $\lambda=\sqrt{-i\omega\tau_{v}}=\alpha a$. The function $\mathcal{Y}_{0}(\omega)$ may be cast in the form 11

where $M_{0}=m_{a}/(6\pi a^{2}\rho)$, and the function $\psi(\lambda)$ tends to zero at large $\lambda$. For a sphere $R=a,\;M_{0}=a/9$, and $\psi(\lambda)=0$. In that case the expression for $\mathcal{Y}_{0}(\omega)$ has just two poles in the complex $\lambda$ plane.

More generally, the function $\psi(\lambda)$ may be written as a continued fraction

A Padé approximant to the exact function is obtained by breaking off at some level. A three-pole approximation to the admittance is found by putting

with suitably chosen coefficients $A_{0},\;A_{1}$. In order to agree with the high-frequency behavior given by Eq. (11) we must have

in the three-pole approximation. Lawrence and Weinbaum 4 have suggested on the basis of their explicit calculations for a spheroid that a good approximation to the admittance $\mathcal{Y}_{0}(\omega)$ is found by use of Eq. (16) with $A_{1}=1$ and the coefficient $A_{0}$ calculated from the Basset coefficient $B$ by use of Eq. (17). In the following we examine this suggestion for a test case, for both an oblate and a prolate spheroid.

The calculation will be based on the solution of the linearized Navier-Stokes equations Eq. (2) after a Fourier transform in time. By axial symmetry the flow velocity and the pressure can be derived from a scalar Stokes stream function. The partial differential equation for the stream function separates in spheroidal coordinates, and the solution can be expressed as an infinite sum of products of radial and angular spheroidal wavefunctions 4 . In the calculation we have used spheroidal wavefunctions as programmed in Mathematica 6.

It has been shown by Lawrence and Weinbaum 16 that the hydrodynamic force $K_{\omega}$ can be evaluated from the behavior of the stream function at large distance from the origin, i.e. from the value of the lowest order multipole moment. The exact set of coupled multipole equations has been derived by Lai and Mockros 3 . It was shown by Lawrence and Weinbaum 4 that the infinite set of linear coupled multipole equations can be solved approximately by truncation and matrix inversion.

It is convenient to define the dimensionless reduced hydrodynamic force $G(\omega)$ by

It follows from Eq. (13) that

We denote the corresponding expression obtained by putting $\psi(\lambda)=0$ in Eq. (14) as

This is a two-pole approximation to the actual function. We call this a two-pole approximation, since the function has two poles in the complex $\sqrt{\omega}$ plane. The deviation from the exact function is characterized by the function $\psi(\lambda)$ in Eq. (14).

We have used the exact multipole equations to calculate the function $\psi(\lambda)$ at a chosen set of positive values $\lambda_{j}$ corresponding to values on the positive imaginary axis of the complex frequency plane. As test cases we have chosen a prolate and an oblate spheroid with aspect ratio $b/a=0.5$. It follows from Table 4-26.1 in Happel and Brenner 15 that the corresponding value of the steady-state friction coefficient is $\zeta(0)=0.9053\times 6\pi\eta a$ in the oblate case, and $\zeta(0)=1.2039\times 6\pi\eta b$ in the prolate case. For values of the aspect ratio $b/a$ closer to unity the correction $\psi(\lambda)$ is very small, and hardly worth calculating. For values of the aspect ratio much less than unity a large number of multipoles must be included to achieve convergence, and the numerical calculation becomes more demanding.

In Table I we list the values $\psi_{j}=\psi(\lambda_{j})$ at the chosen values $\lambda_{1}=0.2,\;\lambda_{2}=0.4,\;\lambda_{3}=0.7,\;\lambda_{4}=1,\;\lambda_{5}=1.5,\;\lambda_{6}=2$, and compare with the corresponding values of the function

corresponding to the approximation of Lawrence and Weinbaum 4 . For the oblate spheroid there is fair agreement, but for the prolate spheroid the LW-approximation yields values with the wrong sign. The numerator in Eq. (21) is negative for the prolate spheroid. It is positive for an oblate spheroid.

In this section we derive approximate expressions for the velocity relaxation function based on the results derived above, and show numerical results for a prolate and an oblate spheroid of aspect ratio $b/a=0.5$. By inversion of the Fourier transform in Eq. (5) we find

The initial value of the velocity relaxation function $\Psi(t)$ is

with effective mass $m^{*}$. The initial value is less than $1/m_{p}$ due to loss of momentum to infinity via longitudinal sound waves 17 . The low frequency behavior of the friction tensor, given by Eq. (9), corresponds to the long-time behavior 18

independent of the shape and size of the spheroid and its mass density $\rho_{p}$.

We write the relaxation function in the form

with dimensionless function $\gamma(t)$ with initial value $\gamma(0)=1$. The function has the one-sided Fourier transform

The simplest estimate of the function $\gamma(t)$ is the exponential $\exp(-t/\tau_{M})$, with mean relaxation time $\tau_{M}=m^{*}/(6\pi\eta R)$, but this violates the long-time behavior given by Eq. (24).

An approximate relaxation function which accounts for the initial value, the mean relaxation time $\tau_{M}$, and the long-time behavior (24), is given by the approximate two-pole expression

where $M=m^{*}/(6\pi a^{2}\rho)$ has the dimension length. In the two-pole approximation the function $\gamma(t)$ takes the approximate form 18

where

and $w(z)=\exp(-z^{2})\mathrm{erfc}(-iz)$.

The two-pole approximation is determined by the effective mass and the friction coefficient at zero frequency. In order to find the actual behavior of the velocity relaxation function, we must take account of the function $\psi(\lambda)$, defined in Eq. (14). The approximation of Lawrence and Weinbaum, corresponding to Eq. (21), leads to a three-pole approximation to the velocity relaxation function, which we denote as $\gamma_{LW}(t)$. More generally, the time-dependence of the reduced velocity relaxation function in an $N$-pole approximation is given by 14

where $\{v_{k}\}$ are the $N$ poles of the one-sided Fourier transform $\Gamma_{N}(\omega)$ in the complex $\lambda$ plane, and $\{R_{k}\}$ are the residues.

The function $\Gamma_{N}(\omega)$ is found numerically from a Padé approximant to the function $\psi(\lambda)$ determined by the values of the reduced hydrodynamic force $G(\omega)$ at selected positive values $\{\lambda_{j}\}$. For example, the six values listed in Table I are sufficient to determine a 5-pole approximation $\gamma_{5}(t)$. In this case the Padé approximant to the function $\psi(\lambda)$ is the ratio of a cubic and a quintic polynomial in $\lambda$.

We find numerically that for aspect ratio $b/a=0.5$ it suffices to include multipoles of order $n=17$, in the notation of Lawrence and Weinbaum 4 . This corresponds to a set of nine coupled multipole equations, since only odd multipoles occur. We have calculated the values of $\psi(\lambda)$ at eight positive values $\{\lambda_{j}\}$, corresponding to a six-pole approximation $\gamma_{6}(t)$. It turns out that in fact two of the residues are negligibly small, so that $\gamma(t)$ is represented accurately by four pole contributions. The larger number of values $\{\lambda_{j}\}$ allow a more accurate determination of the relevant residues and poles. The difference between the approximations $\gamma_{5}(t)$ and $\gamma_{6}(t)$ cannot be seen graphically.

In Fig. 1 we plot the ratio $\gamma_{6}(t)/\gamma_{2p}(t)$ as a function of $\log_{10}(t/\tau_{v})$ for an oblate spheroid with aspect ratio $b/a=0.5$ and vanishing mass density $\rho_{p}=0$, and compare with the approximation of Lawrence and Weinbaum calculated from Eq. (21). In Fig. 2 we present the corresponding plots for a neutrally buoyant spheroid with $\rho_{p}=\rho$. In Figs. 3 and 4 we present similar plots for a prolate spheroid with aspect ratio $b/a=0.5$ and mass densities $\rho_{p}=0$ and $\rho_{p}=\rho$, respectively. It turns out that for the oblate spheroid the Lawrence-Weinbaum approximation agrees qualitatively, but that for the prolate spheroid it is qualitatively wrong. In both cases the deviation from the two-pole approximation $\gamma_{2p}(t)$ is less than one percent.

The frequency-dependent admittance of a prolate or oblate spheroid can be analyzed efficiently by the method of Padé approximants applied to the approximate solution of the exact set of coupled equations for the multipole moments of the hydrodynamic force density. The method leads to an approximate representation of the admittance as a ratio of two polynomials in the square root of frequency. Hence one can evaluate the time-dependent velocity relaxation function as a sum of $w$ functions.

We have applied the method to a prolate or oblate spheroid of aspect ratio $b/a=0.5$ for two values of the spheroid mass density, and have compared the numerical results with those calculated from the approximation proposed by Lawrence and Weinbaum 4 . The comparison is shown in Figs. 1-4. The figures show that for the cases considered the deviation from the two-pole approximation is quite small, and not well described by the approximation of Lawrence and Weinbaum 4 , based on just the Basset coefficient.

The elaborate calculation in terms of spheroidal wavefunctions leads to only a small deviation from the result which may be calculated from the steady state friction coefficient and the effective mass, as given by the classical expressions for a general ellipsoid 2 ^,11 . Presumably an exact calculation for an ellipsoid would lead to the same conclusion in that case. For many practical purposes 19 the two-pole approximation to the admittance, as determined from steady state friction coefficient and effective mass, will be sufficient. The use of only the Basset coefficient for a calculation of the deviations from the two-pole approximation on the intermediate timescale is questionable.

Table of values of the function $\psi(\lambda)$ at selected values of $\lambda$ for an oblate and a prolate spheroid of aspect ratio $b/a=0.5$, as calculated from the set of coupled multipole equations, compared with the Lawrence-Weinbaum approximation given by Eq. (21).