1. Introduction
It is known from potential flow theory that the hydrodynamic force on a solid sphere in arbitrary motion in an incompressible fluid otherwise at rest is equivalent to an increase of the sphere’s mass, an equivalence captured by the label of added, or virtual, mass used for this force (see e.g. Batchelor Reference Batchelor1967; Lighthill Reference Lighthill1986; Landau & Lifshitz Reference Landau and Lifshitz1987). The elementary expression for this force, namely
$\boldsymbol {F}^{am}=-({1}/{2})\rho V \dot {\boldsymbol {w}}$
, in which
$\rho$
is the fluid density,
$V$
the sphere’s volume and
$ \dot {\boldsymbol {w}}$
the sphere’s acceleration, is insufficient for many important applications. An example among many is the Euler–Lagrange computation of disperse multiphase flow in which the particles, far from moving in a quiescent fluid, are typically immersed in a complex flow. It is therefore of some interest to enquire whether formulations of the added mass force of broader applicability exist, even in the limited context of incompressible potential flow. In view of the very different mathematical structure of Stokes vs. potential flow (especially the liner nature of the former vs. the nonlinear nature of the latter), it would be unreasonable to expect a complete generalisation comparable to the one that Faxèn’s theorem provides for the Stokes drag law, but non-trivial potential flows may exist for which such a generalisation may be possible. This is the problem investigated in the present study by considering a fairly general class of potential flows and determining the force acting on a spherical particle in arbitrary motion in them.
To be sure, limited extensions of the added mass force already exist. For example if, in the absence of the sphere, the fluid is subjected to a spatially uniform pressure gradient which puts it into solid-body motion with an ambient velocity
$\boldsymbol {U}^a(t)$
, the previous expression for the force is modified in an obvious way resulting in
${\boldsymbol {F}}^{am}= -({1}/{2})\rho V (\dot {\boldsymbol { w}}-\dot {\boldsymbol {U}}^a)$
(see e.g. Batchelor Reference Batchelor1967). In this case, however, the imposed ambient pressure gradient
$ \boldsymbol{\nabla }p^a$
also exerts a force on the sphere given by
$\boldsymbol {F}^p=-V\boldsymbol{\nabla }p^a$
which, from Euler’s equation of motion, may also be written as
$\boldsymbol {F}^p=\rho V\dot {\boldsymbol {U}}^a$
. Thus, the equation of motion of a sphere of mass
$m$
becomes
\begin{align} m\dot {\boldsymbol {w}} = \boldsymbol {F}^p+\boldsymbol {F}^{am} = \frac {3}{2}\rho V \dot {\boldsymbol {U}}^a-\frac {1}{2}\rho V\dot {\boldsymbol { w}} . \end{align}
Other forces unrelated to the ambient fluid, such as gravity, could readily be included here but, since they would play no role in what follows, they have been omitted for simplicity.
A further extension of the added mass force was given by Taylor (Reference Taylor1928) who showed that, when the ambient fluid velocity
$\boldsymbol {u}^a(\boldsymbol {x},t)$
varies little over spatial scales comparable to the body size, the first term on the right-hand side of (1.1) takes the form
$- V (1 +\alpha ) \boldsymbol{\nabla }p^a$
, in which ‘
$\alpha \rho V$
is the virtual addition to the mass of the body’ (with
$\alpha = {1}/{2}$
for a sphere) and
$p^a$
is the pressure of the ambient flow, namely the flow in the absence of the sphere. Upon using the fluid momentum equation to express
$\boldsymbol{\nabla }p^a$
, (1.1) can then be written as
\begin{align} m\dot {\boldsymbol {w}} = \rho V \frac {{\rm D}\boldsymbol {u}^a}{{\rm D}t} -\frac {1}{2}\rho V\left (\dot {\boldsymbol {w}}-\frac {{\rm D}\boldsymbol { u}^a}{{\rm D}t}\right ) \!, \end{align}
which can be interpreted as showing that, in this case, the added mass force is given by
\begin{align} \boldsymbol {F}^{am} = -\frac {1}{2}\rho V \left (\dot {\boldsymbol {w}}-\frac {{\rm D}\boldsymbol {u}^a}{{\rm D}t}\right )\! . \end{align}
The results of our investigation show that, within the broad class of potential flows that we consider, the previous two extensions are the only ones that can be interpreted in terms of the concept of ‘added mass’. While the fluid force is found to depend linearly on the sphere’s acceleration in the same way shown in (1.2), the no-penetration condition at the sphere’s surface, combined with the quadratic dependence of the pressure on the fluid velocity, modifies the force by introducing a series of terms nonlinear in the spatial derivatives of
$\boldsymbol {u}^a$
. By their very nature, such terms cannot be expressed as contributions to any of the usual time derivatives of the fluid velocity.
Except for the representation used for the potential of the ambient flow and the evaluation of the integral shown in Appendix A, the calculation we present is routine and is summarised in §§ 2–4 which cover the problem set-up, the force due to the ambient flow and the force due to the disturbance flow induced by the particle. Section 5 presents the final result and other pertinent comments, including an example, the calculation of the fluid kinetic energy and a comparison of the main result with Kelvin’s theory of the impulse.
2. Sphere in a general potential flow
If the fluid is inviscid, the particle force that needs to be considered is only due to pressure and is given by
\begin{align} \boldsymbol {F} (\boldsymbol {y}(t),t) = -\oint _S [p]_R\boldsymbol {n}{\rm d}S_x , \end{align}
with
$\boldsymbol {n} =[\boldsymbol {x}-\boldsymbol {y}]_R/R$
the unit normal to the sphere’s surface
$S$
directed into the fluid and
$\boldsymbol {y}(t)$
the instantaneous position of the particle. Here and in the following, the notation
$[\ldots ]_R$
, with
$R$
the radius of the sphere, indicates evaluation on
$S$
and we write
${\rm d}S_x$
as a reminder that integration is over the variable
$\boldsymbol {x}$
implicit in
$p$
and
$\boldsymbol {n}$
while
$\boldsymbol {y}$
is held constant.
By writing the pressure as
$p=p^a+p^d$
, with
$p^a$
the pressure of the ambient flow in the absence of the particle and
$p^d$
the pressure disturbance caused by the particle, we can decompose the force into an ambient force
$ \boldsymbol {F}^a$
and a disturbance force
$\boldsymbol {F}^d$
as
\begin{align} \boldsymbol {F} = \boldsymbol {F}^a + \boldsymbol {F}^d \, , \qquad \textrm {with} \qquad \boldsymbol {F}^{a,d} = -\oint _S \big[p^{a,d} \big]_R\boldsymbol {n}{\rm d}S_x . \end{align}
In an irrotational flow the pressure can be expressed in terms of the Bernoulli integral
\begin{align} \partial _t\phi +\frac {1}{2} u^2 +\frac {p}{\rho } = f(t) , \end{align}
in which
$\partial _t\equiv \partial /\partial t$
denotes the partial time derivative and
$f(t)$
is a function of time which need not concern us as its contribution to the force vanishes. Upon writing
\begin{align} -\frac { p^a}{\rho } =\partial _t\phi ^a +\frac {1}{2} (u^a)^2-f(t) \,,\qquad -\frac { p^d}{\rho } =\partial _t(\phi -\phi ^a) +\frac {1}{2} \big(u^2-(u^a)^2 \big) , \end{align}
with
$\phi ^a$
the potential of the undisturbed ambient flow and
$\boldsymbol {u}^a=\boldsymbol{\nabla} _{\boldsymbol {x}}\phi ^a(\boldsymbol { x},t)$
the corresponding velocity, and substituting into (2.2), we can express the two components of the force in terms of the velocity potentials of the ambient and disturbed flows.
By writing
$\phi ^a(\boldsymbol {x},t)=\phi ^a(\boldsymbol {y}+(\boldsymbol {x}-\boldsymbol {y}),t)$
and expanding in Taylor series around
$\boldsymbol {y}$
, the ambient potential can be conveniently expressed as
\begin{align} \phi ^a(\boldsymbol {x},t) = \sum _{l=0}^\infty {1\over l!} \varPhi _l , \end{align}
in which
$\varPhi _l$
is defined by
\begin{align} \varPhi _l= [(\boldsymbol {x}-\boldsymbol {y})\boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l)} \phi ^a(\boldsymbol {y},t)\, \equiv \, (x-y)_i(x-y)_j\ldots (x-y)_k\partial _i\partial _j\ldots \partial _k \phi ^a(\boldsymbol {y},t), \end{align}
with
$\partial _i\equiv \partial / \partial y_i$
and the derivatives evaluated at y; for the special case
$l=0$
,
$\varPhi _0=\phi ^a(\boldsymbol {y},t)$
. Here and in the following, unsubscripted derivatives are with respect to the position
$\boldsymbol {y}$
of the sphere’s centre; derivatives with respect to the field point
$\boldsymbol {x}$
are indicated by
$\boldsymbol{\nabla} _{\boldsymbol {x}}$
.
Use of the Taylor-series representation (2.5) of the ambient potential raises obvious convergence questions. In most of this paper this representation is used only in the neighbourhood of the sphere and therefore, while somewhat limiting in the generality of the conclusions that follow, this assumed form for
$\phi ^a$
is not overly restrictive in the sense that it still allows for non-trivial potential flows, most obviously when the series terminates. When it does not, it may be expected that its radius of convergence would decrease with increasing
$|\boldsymbol {x}-\boldsymbol {y}|$
. This may present a problem for the calculation of the excess kinetic energy of the flow in § 5 which, however, would still be justified when the series terminates and possibly in other situations as well.
The ambient velocity corresponding to the potential (2.5) is
\begin{align} \boldsymbol {u}^a(\boldsymbol {x},t)= \boldsymbol{\nabla} _{\boldsymbol {x}}\phi ^a=\sum _{l=0}^\infty \frac {1}{l!} \boldsymbol{\nabla} _{\boldsymbol {x}}\varPhi _l = \sum _{l=1}^\infty \frac {1}{(l-1)!} [(\boldsymbol {x}-\boldsymbol {y})\boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l-1)} \boldsymbol{\nabla }\phi ^a(\boldsymbol {y},t) , \end{align}
and, on the particle’s surface, it reduces to
\begin{align} [\boldsymbol {u}^a]_R = \sum _{l=0}^\infty \frac {R^l}{l!} \left [\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ]^{(l)} \boldsymbol{\nabla }\phi ^a(\boldsymbol {y},t) \, . \end{align}
Here,
$ [\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l)} = n_i n_ j \ldots n_k \partial _i\partial _j\ldots \partial _k$
is the tensor product of
$l$
normals with
$l$
derivatives.
We write the disturbance potential
$\phi ^d$
in the form
\begin{align} \phi ^d(\boldsymbol {x},t)=\phi -\phi ^a= {\mathcal W} + \psi ^d , \end{align}
with the first term representing the potential of a sphere moving in a fluid at rest, namely
\begin{align} {\mathcal W} = -{1\over 2} \frac {R^3}{|\boldsymbol {x}-\boldsymbol {y}|^3}\,\boldsymbol {w}\boldsymbol{\cdot }\left (\boldsymbol {x}-\boldsymbol {y}\right ) \,,\qquad \textrm {and} \qquad \psi ^d = \sum _{l=0}^\infty \frac {l}{(l+1)!} \left (\frac {R}{|\boldsymbol {x}-\boldsymbol {y}|}\right )^{2l+1} \varPhi _l . \end{align}
The corresponding disturbance velocity on the sphere’s surface is
\begin{eqnarray} \big[\boldsymbol {u}^d \big]_R &=& [\boldsymbol{\nabla} _{\boldsymbol {x}} {\mathcal W}]_R + \big[\boldsymbol{\nabla} _{\boldsymbol {x}} \psi ^d \big]_R= \frac {3}{2} (\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol { n})\boldsymbol {n}-\frac {1}{2}\boldsymbol {w}\nonumber \\ && \quad +\sum _{l=1}^\infty \frac {l}{(l+1)!} R^{l-1}\big (l[\boldsymbol {n}\boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l-1)}\boldsymbol{\nabla }\phi ^a-(2l+1)\boldsymbol {n}[\boldsymbol {n}\boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l)}\phi ^a\big ). \end{eqnarray}
It is readily shown that, on the surface of the particle,
\begin{align} \boldsymbol {n}\boldsymbol{\cdot } \big[\boldsymbol{\nabla }\big(\phi ^a+\psi ^d \big) \big]_R=0 \,, \qquad \textrm {so that} \qquad \big[\boldsymbol {n}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol {x}}\phi ^d \big]_R=[\boldsymbol {n}\boldsymbol{\cdot }\boldsymbol{\nabla }{\mathcal W}]_R= \boldsymbol {n}\boldsymbol{\cdot }\boldsymbol {w} \, . \end{align}
Furthermore, it is evident that
$\phi ^d \rightarrow 0$
far from the particle.
3. Ambient force
Upon using the Bernoulli integral, the ambient force is
\begin{align} \boldsymbol {F}^a = \rho \oint _S \left (\partial _t\phi ^a+\frac {1}{2}(u^a)^2\right )\boldsymbol {n}\,{\rm d}S_x \, . \end{align}
The time derivative of the ambient potential
$\phi ^a$
is given by (2.5) with
$\varPhi _l$
replaced by
\begin{align} \partial _t\varPhi _l = -l [(\boldsymbol {x}-\boldsymbol {y})\boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l-1)}\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol{\nabla }\phi ^a + [(\boldsymbol {x}-\boldsymbol {y})\boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l)} (\partial _t \phi ^a+\boldsymbol {w} \boldsymbol{\cdot }\boldsymbol{\nabla }\phi ^a) , \end{align}
from which, on the sphere’s surface, the time derivative of the ambient potential follows as
\begin{align} [\partial _t \phi ^a]_R= \sum _{l=0}^\infty \frac {1}{l!} R^{l-1}\big [-l [\boldsymbol {n}\boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l-1)}\boldsymbol { w}\boldsymbol{\cdot }\boldsymbol{\nabla }\phi ^a +R[\boldsymbol {n}\boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l)} (\partial _t \phi ^a+\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol{\nabla }\phi ^a) \big ]. \end{align}
The contribution to the integral (3.1) over the sphere’s surface is readily found to be
\begin{align} \oint _S [\partial _t\phi ^a]_R \boldsymbol {n}{\rm d}S_x = V [\partial _t \boldsymbol {u}^a+(\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol {u}^a] - V (\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol {u}^a =V\partial _t \boldsymbol {u}^a(\boldsymbol {y},t). \end{align}
This simple result follows from the fact that integrals over a sphere of an odd number of normals vanish. Since the integrand explicitly shows one normal, only terms of
$ [\partial _t\phi ^a]_R$
with an odd number of normals can contribute. However, of these, all terms with more than a single normal must also vanish as they will necessarily contain factors proportional to
${\nabla} ^2\phi ^a = 0$
; further comments on this point will be found in Appendix A.
With
$d_l=R^l/l!$
, the quadratic contribution to the integral in (3.1) is
\begin{align} \frac {1}{2} \oint _S (u^a)^2n_i {\rm d}S_x = \frac {1}{2}\sum _{l=0}^\infty d_l \sum _{k=0}^\infty d_k\oint _S \big ( \left [\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ]^{(l)} \boldsymbol{\nabla }\phi ^a\big ) \boldsymbol{\cdot }\big ( \left [\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ]^{(k)} \boldsymbol{\nabla }\phi ^a\big ) n_i{\rm d}S_x . \end{align}
For a given
$l$
, only terms with
$k = l\pm 1$
give a non-vanishing contribution to the integral as all others include
${\nabla} ^2\phi ^a$
or its derivatives (see Appendix A). Furthermore, since the integrals are symmetric in the indices
$l$
and
$k$
, for any value
$l^*$
of
$l$
, the contribution of the term
$l=l^*$
and
$k =l^*+1$
is equal to the contribution of the term
$l=l^*+1$
and
$k =l^*$
. Thus the double sum reduces to a single sum in which all terms appear twice
\begin{align} \frac {1}{2} \oint _S (u^a)^2n_i {\rm d}S_x = \frac {1}{2} \sum _{l=0}^\infty 2d_ld_{l+1} \oint _S \left ( \left [\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ]^{(l)} \boldsymbol{\nabla }\phi ^a\right ) \boldsymbol{\cdot }\left ( \left [\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ]^{(l+1)} \boldsymbol{\nabla }\phi ^a\right ) n_i{\rm d}S_x . \end{align}
Integrals of this form are evaluated in Appendix A and, upon using the result given in (A8), we find an expression that can be written as
\begin{align} \frac {1}{2} \oint _S (u^a)^2\boldsymbol {n} {\rm d}S_x = V (\boldsymbol {u}^a\boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u}^a (\boldsymbol {y},t) +3V \sum _{l=1}^\infty \frac {R^{2l}}{l!\,(2l+3)!!} \big\{(\partial ^l \boldsymbol {u}^a)\boldsymbol{\cdot }\boldsymbol{\nabla } \big\} \big(\partial ^l \boldsymbol { u}^a \big) , \end{align}
in which
$(2l+3)!!=(2l+3)(2l+1) \ldots 3\boldsymbol{\cdot }1$
, the double factorial of
$2l+3$
, is the product of all odd numbers less than, or equal to,
$2l+3$
. In (3.7) we have introduced the notation
\begin{align} \big [\big\{\big(\partial ^l \boldsymbol {u}^a \big)\boldsymbol{\cdot }\boldsymbol{\nabla }\big\} \big(\partial ^l \boldsymbol {u}^a \big)\big ]_i = (\underbrace {\partial _a\partial _b \ldots \partial _c}_{l \,\, {terms}} u_{\kern-1.5pt j}^a)(\partial _j \underbrace {\partial _a\partial _b \ldots \partial _c}_{l \,\,{terms}}u_i ^a), \end{align}
to be understood in the sense of a normal tensor product. This abbreviated writing is reminiscent of that used for the nonlinear term of the convective derivative, namely
$(u_{\kern-1.5pt j}^a\partial _j)u_i^a$
, with the difference that the velocity
$u_i^a$
is replaced by its multi-gradient
$\partial ^lu^a_i=\partial _a\partial _b \ldots \partial _c u_i ^a$
. With this notation, for
$l=0$
, the familiar nonlinear contribution to the convective derivative appearing in the first term on the right-hand side of (3.7) is recovered.
Upon noting that, for potential flow,
$\partial _j u^a_i = \partial _i u_{\kern-1.5pt j}^a$
, this result can also be written as
\begin{align} \frac {1}{2}\oint _S (u^a)^2\boldsymbol {n} {\rm d}S_x = V \left [ (\boldsymbol {u}^a\boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u}^a +\boldsymbol{\nabla }\varPi ^a \right ], \end{align}
in which the potential per unit displaced fluid mass
$\varPi ^a$
is given by
\begin{align} \varPi ^a=\frac {3}{2} \sum _{l=1}^\infty \frac {R^{2l}}{l!\,(2l+3)!!} \big(\partial ^l u_{\kern-1.5pt j}^a \big)^2 =\frac {R^2}{10}(\boldsymbol{\nabla }\boldsymbol {u}^a):(\boldsymbol{\nabla }\boldsymbol {u}^a) +\frac {R^4}{140} (\boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol { u}^a)\stackrel {.}{:}(\boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol {u}^a)+\ldots . \end{align}
Here,
$(\partial ^l u_{\kern-1.5pt j}^a)^2=(\partial _a\partial _b \ldots \partial _c u_{\kern-1.5pt j}^a)\,(\partial _a\partial _b \ldots \partial _c u_{\kern-1.5pt j}^a)$
in the standard sense of tensor products. With (3.4) and (3.9) the ambient force is then found in the form
\begin{align} \boldsymbol {F}^a = \rho V \big[\partial _t \boldsymbol {u}^a+(\boldsymbol {u}^a\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol {u}^a \big] +\rho V \boldsymbol{\nabla }\varPi ^a = \rho V\left (\frac {{\rm D}\boldsymbol {u}^a}{{\rm D}t} +\boldsymbol{\nabla }\varPi ^a\right ). \qquad \end{align}
This result is a generalisation of the pressure-gradient force shown in the first term on the right-hand side of (1.2). It may be of some interest to note that
$\boldsymbol {F}^a$
can also be written as
$\rho V\boldsymbol{\nabla }(\partial _t\phi ^a+({1}/{2}) \boldsymbol {u}^a\boldsymbol{\cdot }\boldsymbol {u}^a +\varPi ^a)$
. Not surprisingly, the hydrodynamic force acting on the particle in a general inviscid flow is found to be conservative.
The infinite series defining
$\varPi ^a$
shows that all the scales of the ambient flow contribute to the force, as expected. With the very rough estimate
$ R^{2l}(\partial ^l u_{\kern-1.5pt j}^a)^2 \sim (R/L)^{2l}U^2$
, with
$U$
and
$L$
representative values of the ambient velocity and the related spatial scale, for
$R\lt L$
we see that the effect of
$\varPi ^a$
is small and decreases with increasing
$l$
. Thus we recover Taylor’s result. For
$R/L \gt 5$
or so, we can use Stirling’s formula to obtain a rough approximation to the integer
$l_{max}$
which gives the largest contribution to the sum. Unsurprisingly, the result is
$l_{max}\sim R/L$
. The coefficients of the summation decrease very rapidly as
$l$
increases beyond a value of this order indicating that the finer structure of the ambient flow is inconsequential for the ambient force.
4. Disturbance force
On the sphere’s surface the time derivative of the disturbance potential (2.9) is the sum of
\begin{align} [\partial _t {\mathcal W}]_R= \frac {1}{2}\big (w^2-R\dot {\boldsymbol {w}}\boldsymbol{\cdot }\boldsymbol {n}\big )-\frac {3}{2}(\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol {n})^2 , \end{align}
and
\begin{align} [\partial _t\psi ^d]_R & = \sum _{l=0}^\infty \frac {l R^{l-1}}{(l+1)!} \big ((2l+1)(\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol { n})[\boldsymbol {n}\boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l)} \phi ^a - l [\boldsymbol {n}\boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l-1)}(\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol{\nabla })\phi ^a \big . \nonumber \\& \qquad \big. +R[\boldsymbol {n}\boldsymbol{\cdot }\boldsymbol{\nabla }]^{(l)}(\partial _t\phi ^a+\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol{\nabla }\phi ^a)\big ) . \end{align}
The contribution of these two terms to the integral analogous to that shown in (3.4) is
\begin{align} \oint _S \big[\partial _t \phi ^d \big]\boldsymbol {n}{\rm d}S_x = \frac {1}{2} V \big[\partial _t \boldsymbol {u}^a +(\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol { u}^a-\dot {\boldsymbol {w}}\big ] . \end{align}
To calculate the nonlinear contribution it is convenient to note that, in view of the vanishing of the normal component of
$\boldsymbol{\nabla} _{\boldsymbol {x}} (\phi ^a+\psi ^d)$
on the particle’s surface, the square of the disturbance velocity may be written as
\begin{align} \big[(u^d)^2 \big]_R & = \big[(\boldsymbol{\nabla} _{\boldsymbol {x}}{\mathcal W})^2 \big]_R + 2 \big[ \big(\boldsymbol{\nabla} _{\boldsymbol {x}} \phi ^a+\boldsymbol{\nabla} _{\boldsymbol {x}} \psi ^d \big)_{tg}\boldsymbol{\cdot }(\boldsymbol{\nabla} _{\boldsymbol {x}} {\mathcal W})_{tg} \big]_R \nonumber \\& \quad + \big[\big(\boldsymbol{\nabla} _{\boldsymbol {x}} \phi ^a+\boldsymbol{\nabla} _{\boldsymbol {x}} \psi ^d \big)_{tg}^2 \big]_R - \big[(u^a)^2 \big]_R, \end{align}
with
$(\ldots )_{tg}$
indicating the component tangent to the sphere’s surface. The contribution of this term to the force is to be found by performing an integration similar to that shown on the left-hand side of (3.5) with
$(u^a)^2$
replaced by
$(u^d)^2$
. Once multiplied by
$n_i$
, the integral of the first term of (4.4) includes an odd number of normals and, therefore, vanishes. The integral of the last term has already been calculated in (3.7). Use of the identity
$[\boldsymbol{\nabla} _{\boldsymbol {x}} (\phi ^a+\psi ^d)]_{tg}\boldsymbol{\cdot }(\boldsymbol{\nabla} _{\boldsymbol {x}} {\mathcal W})_{tg} = [\boldsymbol { n}\times \boldsymbol{\nabla} _{\boldsymbol {x}} (\phi ^a+\psi ^d)]\boldsymbol{\cdot }[\boldsymbol {n}\times \boldsymbol{\nabla} _{\boldsymbol {x}} {\mathcal W}]$
and similarly for the third term simplifies the remaining integrations. The contribution of the second term of (4.4) is then found to be
\begin{align} \frac {1}{2}\oint _S \big[2\big(\boldsymbol{\nabla} _{\boldsymbol {x}} \phi ^a+\boldsymbol{\nabla} _{\boldsymbol {x}} \psi ^d \big)_{tg}\boldsymbol{\cdot }(\boldsymbol{\nabla} _{\boldsymbol {x}} {\mathcal W})_{tg} \big]_R n_i {\rm d}S_x =-\frac {1}{2}V(\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol{\nabla }) u_i^a(\boldsymbol {y},t) \, . \end{align}
With the notation introduced before in (3.7), the contribution of the third term of (4.4) is
\begin{align} \frac {1}{2}\oint _S \big[\big(\boldsymbol{\nabla} _{\boldsymbol {x}} \phi ^a+\boldsymbol{\nabla} _{\boldsymbol {x}} \psi ^d \big)_{tg}^2 \big]_R n_i {\rm d}S_x & = \frac {3}{2} V(\boldsymbol {u}^a\boldsymbol{\cdot }\boldsymbol{\nabla })u_i^a \nonumber \\ & \quad + 3V\sum _{l=1}^\infty \frac {R^{2l}}{(l+2) l!(2l+1)!!} \big\{(\partial ^l \boldsymbol { u}^a)\boldsymbol{\cdot }\boldsymbol{\nabla } \big\} \big(\partial ^l u_i^a \big) . \end{align}
After subtraction of (3.7) as per the last term of (4.4), the combination of (4.5) and (4.6) gives the contribution of the square of the disturbance velocity to the force as
\begin{align} \frac {1}{2} \oint _S \big(u^d \big)^2\boldsymbol {n} {\rm d}S_x &= -\frac {1}{2}V \big[(\boldsymbol {w}-\boldsymbol {u}^a)\boldsymbol{\cdot }\boldsymbol{\nabla }\big] \boldsymbol { u}^a \nonumber \\ & \quad +3V \sum _{l=1}^\infty \frac {l+1}{l+2} \frac {R^{2l}}{l!(2l+3)!!} \big\{\big(\partial ^l \boldsymbol {u}^a \big)\boldsymbol{\cdot }\boldsymbol{\nabla }\big\}\big(\partial ^l \boldsymbol {u}^a\big). \end{align}
As before, the summation may be written as
$V\boldsymbol{\nabla }\varPi ^d$
with
$\varPi ^d$
given by
\begin{align} \varPi ^d = \frac {3}{2}\sum _{l=1}^\infty \frac {(l+1)R^{2l}}{(l+2) l!(2l+3)!!} \big(\partial ^l u_{\kern-1.5pt j}^a \big)^2 = \varPi ^a-\frac {3}{2}\sum _{l=1}^\infty \frac {R^{2l}}{(l+2) l!(2l+3)!!} \big(\partial ^l u_{\kern-1.5pt j}^a \big)^2 \, . \end{align}
Upon combining (4.3) and (4.7), the disturbance force is finally found to be
\begin{align} \boldsymbol {F}^d = -\frac {1}{2} \rho V\left (\dot {\boldsymbol {w}}- \frac {{\rm D}\boldsymbol {u}^a}{{\rm D}t}\right ) +\rho V\boldsymbol{\nabla }\kern-1pt\varPi ^d . \end{align}
We see here the generalisation of (1.3) to a general potential flow.
The coefficients of the potential
$\varPi ^d$
differ from those of
$\varPi ^a$
only by the factor
$(l+1)/(l+2)$
which is close to 1 when the previous estimate
$l_{max} \sim R/L$
is applicable. Thus, the comments made at the end of § 3 with respect to
$l_{max}$
apply to
$\varPi ^d$
as well.
5. The total fluid force and concluding comments
Upon adding (3.11) and (4.9) we find the following extension of the fluid force on the particle given in (1.2) to the case of a general potential flow in the form:
\begin{align} \boldsymbol {F}= \rho V \frac {{\rm D}\boldsymbol {u}^a}{{\rm D}t} -\frac {1}{2} \rho V\left (\dot {\boldsymbol {w}}- \frac {{\rm D}\boldsymbol { u}^a}{{\rm D}t}\right ) +\rho\kern-0.5pt V\kern-0.5pt\boldsymbol{\nabla }\kern-0.5pt\varPi , \end{align}
in which
\begin{align} \varPi &= \varPi ^a+\varPi ^d = \frac {3}{2} \sum _{l=1}^\infty \frac {l+1}{(l+2)!\,(2l+1)!!}R^{2l} \big(\partial ^l u_{\kern-1.5pt j}^a \big)^2 \nonumber \\&=\frac {R^2}{6} (\boldsymbol{\nabla }\boldsymbol {u}^a):(\boldsymbol{\nabla }\boldsymbol {u}^a) +\frac {R^4}{80} (\boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol {u}^a) \stackrel {\boldsymbol{\cdot }}{:} (\boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol {u}^a) + \frac {R^6}{2100} (\boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol {u}^a) \stackrel {:}{:} (\boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol {u}^a) +\ldots . \end{align}
This explicit expression of the first few terms of
$\varPi$
may be used to have an idea of the magnitude of the error affecting the conventional added mass force in specific cases.
It is interesting to note that, if present, the term corresponding to
$l=0$
in the summation defining
$\varPi$
would be
$( {3}/{4}) \boldsymbol {u}^a\boldsymbol{\cdot }\boldsymbol {u}^a$
with a gradient equal to
$( {3}/{2}) (\boldsymbol {u}^a \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol {u}^a$
, namely the sum of the two nonlinear terms implicit in
${\rm D}\boldsymbol {u}^a/{\rm D}t$
in (5.1). On the basis of this remark, when the contribution of the summation is negligible, the form (1.2) might look like an accident although, of course, the ‘accident’ is dictated by two important physical arguments. First, when
$R$
is much smaller than the scale of variation of
$\phi ^a$
, the fluid in the neighbourhood of the particle is essentially in uniform motion with velocity
$\boldsymbol {U}^a(t)\simeq \boldsymbol {u}^a(\boldsymbol {y}(t),t)$
and, by the principle of relativity, the simple expression
$\boldsymbol {F}^{am}=-({1}/{2})\rho V \dot {\boldsymbol {w}}$
for the ambient fluid at rest must become
$\boldsymbol {F}^{am}=-({1}/{2})\rho V (\dot {\boldsymbol {w}}-\dot {\boldsymbol { U}}^a)$
, augmented by the effect of the apparent force
$\rho V \dot {\boldsymbol {U}}^a$
arising from the use of a non-inertial frame of reference. Secondly, neither
$\partial _t\boldsymbol {u}^a$
nor
$(\boldsymbol {u}^a\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol {u}^a$
are Galilean invariants, while their combination into
${\rm D}\boldsymbol {u}^a/{\rm D}t$
is. The term corresponding to
$l=0$
is therefore, in a sense, special and its consideration jointly with
$\partial _t\boldsymbol {u}^a$
is very well justified and, indeed, necessary. This consideration is not applicable to
$\varPi$
as all the terms of the summation include gradients of the velocity and are therefore already Galilean invariant. This property prevents the possibility of seeing
$\varPi$
as part of an added mass force, a concept that the previous result makes inapplicable in the presence of velocity scales comparable to the sphere’s radius.
Some greater insight into the significance of our result can be gained with an explicit example. We consider a sphere of radius
$R$
at rest in the flow generated by a simple point source at the origin with potential
$\phi ^a=-C/r$
;
$C$
is a constant and
$r$
is the distance from the origin. A particularly favourable feature of this example is that, in this case, the series (5.2) can be summed in closed form with the result
\begin{align} \varPi (r)=\frac {3C^2}{2R^4}\left [ \frac {R^2}{r^2} \frac {2r^4-R^2r^2+R^4}{2r^2(r^2-R^2)} +\ln \left(1-\frac {R^2}{r^2} \right)\right ] , \end{align}
in which
$r$
is the distance of the sphere’s centre from the source. In ascending order, figure 1 shows, as functions of
$r/R$
, graphs of
$\rho\kern-0.5pt V\kern-0.5pt \partial _r\kern-0.5pt\varPi$
(lowest blue line), the conventional added mass force
$( {1}/{2})\rho V {\rm D} u^a/{\rm D}t$
(red), the pressure-gradient force
$\rho V {\rm D} u^a/{\rm D}t$
(black) and the entire Bernoulli force (5.1) as obtained by adding up the terms represented by the other lines (light blue). The quantities plotted are the radial components of each term and each one is normalised by division by
$\rho V {\rm D} u^a/{\rm D}t$
with
$u^a$
the radial velocity induced by the source. The crosses on the top light-blue line are numerical results obtained by numerical integration over the sphere’s surface of the fluid pressure as obtained from the Bernoulli theorem. The excellent agreement between the numerical results and the expression (5.1) supports the accuracy of the latter. The force correction (5.3) is singular for
$r=R$
and decays to zero as
$r$
increases. The decay is, however, fairly slow. The simple result
$({3}/{2})\rho V {\rm D} u^a/{\rm D}t$
incurs an error of nearly 45 % for
$r/R=2$
, of approximately 21 % for
$r/R=3$
and of approximately 12 % for
$r/R=4$
.
The components of the pressure force on a sphere of radius
$R$
immersed in a simple source flow as a function of the distance
$r$
of the sphere’s centre from the source; all components are normalised by division by
$\rho V {\rm D} u^a/{\rm D}t$
, with
$u^a$
the radial velocity component of the source flow. In ascending order, the lowest line (blue) shows
$\rho V \partial _r\varPi$
, the next line
$({1}/{2})\rho V {\rm D} u^a/{\rm D}t$
(red) the one above that
$\rho V {\rm D} u^a/{\rm D}t$
(black) and the top line the entire force (5.1). The marks on this line are the result of a numerical evaluation of the pressure force on the sphere.
It is of course legitimate to ask what the real value of our result (5.1) is when it comes to the practical but very important problem of modelling the force on a particle not in potential, but in real viscous flow. The current practice in Euler–Lagrange simulations of disperse flows is to combine the conventional added mass force from potential flow (1.2) with drag and, possibly, the memory force. It is thought that, in this way, a reasonable approximation to the fluid force acting on a suspended particle can be obtained. As long as the flow does not separate and it varies slowly over the particle’s radius, the conventional added mass force may be expected to give a good representation of the inertial component of the fluid force and, indeed, this expectation is supported by numerical work (see e.g. Howe Reference Howe1995; Magnaudet & Eames Reference Magnaudet and Eames2000; Mougin & Magnaudet Reference Mougin and Magnaudet2002; Bagchi & Balachandar Reference Bagchi and Balachandar2003). We may expect our contribution to be similarly useful for an unseparated flow, with the advantage of extending the known result to flows varying on an arbitrarily small spatial scale. Separation, however, changes the pressure distribution behind the particle in a major way and the above commonly adopted procedure becomes questionable. While our result is not subjected to scale limitations, it is still based on unseparated potential flow and, therefore, it can be fairly criticised on the same basis. Rather than in an ability to counter these criticisms, we believe that the real value of our result is a fundamental clarification of the nature of the added mass force.
If the term
$({3}/{4}) \boldsymbol {u}^a\boldsymbol{\cdot }\boldsymbol {u}^a$
is included in the summation defining the potential by writing it as
$\tilde {\varPi }$
to signify a summation starting from
$l=0$
rather than
$l=1$
, (5.1) can be written in an equivalent way which fully exhibits the potential nature of the fluid force on the particle, namely
\begin{align} \boldsymbol {F}=\rho V \boldsymbol{\nabla }\left ( \frac {3}{2}\partial _t\phi ^a + \tilde {\varPi }\right )-\frac {1}{2} \rho V\dot {\boldsymbol { w}} . \end{align}
The force on the sphere can also be viewed from the point of view of Kelvin’s impulse, defined by (see e.g. Lighthill Reference Lighthill1986)
\begin{align} \boldsymbol {J} = -\oint _S \phi \,\boldsymbol {n} {\rm d}S_x \, . \end{align}
Upon using (2.5) and (2.10), this quantity is readily found to be given by
$\boldsymbol { J}=({1}/{2})\rho V \boldsymbol {w}(t)-({3}/{2})\rho V {\boldsymbol {u}}^a(\boldsymbol {y}(t),t)$
. It can be shown that, in terms of
$\boldsymbol {J}$
, the equation of motion for the particle can be written as (see e.g. Blake & Cerone Reference Blake and Cerone1982)
\begin{align} \frac {d}{{\rm d}t}\left (m\boldsymbol {w}+\boldsymbol {J}\right ) = \rho \oint _S \left (\frac {1}{2}u^2\boldsymbol {n}- (\boldsymbol {n}\boldsymbol{\cdot }\boldsymbol {w})\boldsymbol {u}\right ) {\rm d}S_x. \end{align}
The first term on the right-hand side can be found by adding (3.7), (4.5) and (4.7) and the second term is readily calculated to be
$V(\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol {u}^a$
. The impulse equation (5.6) then becomes
\begin{align} m{\dot {\boldsymbol {w}}} & +\left ( \frac {1}{2}\rho V \dot {\boldsymbol {w}} -\frac {3}{2}\rho V \big[\partial _{t}{\boldsymbol {u}}^a+(\boldsymbol {w}\boldsymbol{\cdot }\boldsymbol{\nabla }){\boldsymbol {u}}^a \big] \right ) = \frac {3}{2}\rho V[({\boldsymbol {u}}^a-\boldsymbol {w})\boldsymbol{\cdot }\boldsymbol{\nabla }] \boldsymbol {u}^a \nonumber \\ & +\rho V \boldsymbol{\nabla }\varPi . \end{align}
The expression for the hydrodynamic force implicit in this result is equivalent to (5.1). If the particle moves in a fluid at rest (or in uniform motion), all the terms of this equation except the first two (or, respectively, three) vanish, so that
$-{\rm d}\boldsymbol {J}/{\rm d}t$
can be identified with the hydrodynamic force on the particle and, therefore,
$+{\rm d}\boldsymbol {J}/{\rm d}t$
with the rate of change of the fluid momentum as noted by Kelvin. When the surrounding fluid motion is also influenced by processes unrelated to the particle, however, this simple interpretation of the impulse is no longer possible as indicated by the appearance, on the left-hand side of (5.7), of the convective derivative with the particle’s velocity
$\boldsymbol {w}$
rather than the fluid’s velocity
${\boldsymbol {u}}^a$
. In this more general situation
$\boldsymbol {J}$
, which only depends on the potential as seen by the particle, is not sufficient to determine the fluid force and the terms on the right-hand side of (5.7) become necessary to restore the correct physics.
When the ambient flow is unbounded in all directions, our explicit expression for the potential permits us to obtain the fluid kinetic energy in excess of that in the absence of the sphere by calculating
\begin{align} {\mathcal K}^d = \frac {1}{2}\rho \int _{V_\infty -V} \big\{\big[(\boldsymbol{\nabla }\big(\phi ^a+ \phi ^d \big) \big]^2-(\boldsymbol{\nabla }\phi ^a)^2 \big\} {\rm d}V - \frac {1}{2}\rho \int _V (\boldsymbol{\nabla }\phi ^a)^2 {\rm d}V, \end{align}
in which
$V_\infty$
is a large sphere, concentric with the particle, the radius
$R_\infty$
of which is eventually allowed to become infinite. (The last term needs to be subtracted because the undisturbed kinetic energy is found by integrating over the entire
$V_\infty$
, not just
$V_\infty -V$
.) Use of the previous expressions for the potentials in the calculation of the first integral presupposes that the radius of convergence of the Taylor series used to represent the ambient flow in (2.5) is infinite and that the ambient flow has no singularities at finite distance from the sphere. If these hypotheses are satisfied, it is found that the coefficients of terms involving
$R_\infty$
vanish identically so that the result is finite, and it is found to be
\begin{align} {\mathcal K}^d & = \frac {1}{4}\rho V(\boldsymbol {w}-\boldsymbol {u}^a)^2 - \frac {1}{2}\rho V \boldsymbol {u}^a\boldsymbol{\cdot }\boldsymbol {u}^a -3\rho V\sum _{l=1}^\infty \frac {l+1}{(l+2)!\,(2l+3)!!}R^{2l} \big(\partial ^l u_{\kern-1.5pt j}^a \big)^2 \nonumber \\& = \frac {1}{4}\rho V(\boldsymbol {w}-\boldsymbol {u}^a)^2 - \frac {1}{2}\rho V \boldsymbol {u}^a\boldsymbol{\cdot }\boldsymbol {u}^a \nonumber \\ &\quad -\rho V\left ( \frac {R^2}{15} (\boldsymbol{\nabla }\boldsymbol {u}^a):(\boldsymbol{\nabla }\boldsymbol {u}^a)+\frac {R^4}{280} (\boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol { u}^a)\stackrel {.}{:}(\boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol {u}^a) +\ldots \right ). \end{align}
If the ambient flow were uniform with velocity
$\boldsymbol {u}^a(\boldsymbol {y},t)=\boldsymbol {U}^a(t)$
the summation would be identically zero. The first term would then be the kinetic energy of the flow induced by the sphere’s motion relative to the surrounding ambient flow and the second term would correspond to the reduction of the kinetic energy due to the volume occupied by the sphere as indicated by the last term of (5.8). The summation accounts for the non-uniformity of the ambient flow and the nonlinear interactions among the various velocity scales. The summation in (5.9) differs from that defining
$\varPi$
in (5.2) by an extra factor
$(2l+3)$
in the denominator of the coefficients, which promotes a faster convergence of the series and diminishes the contribution of the smaller scales.
Acknowledgements
The authors are very grateful to Professor L.v. Wijngaarden (University of Twente, The Netherlands) for several very useful comments.
Funding
This work was supported in part by the U.S. Department of Energy, National Nuclear Security Administration, Advanced Simulation and Computing Program, as a Cooperative Agreement under the Predictive Science Academic Alliance Program, under Contract No. DE-NA0004263.
Declaration of interests
The authors report no conflicts of interest.
Appendix A
In this appendix we show the explicit evaluation of integrals such as that appearing in (3.5) and related matters. For clarity let us rewrite the integral in question as
\begin{eqnarray} I_i^{(l,k)} = I_i^{(k,l)} = \oint (\underbrace {n_an_b\ldots n_c}_{l \,\,{terms}} ) (\underbrace {n_m n_p\ldots n_q}_{k \,\,{terms}}) n_i{\rm d}S_x \big(\partial ^l\phi ^a \big) \big(\partial ^k \phi ^a \big) , \end{eqnarray}
with
$\partial ^l\phi ^a= \partial _a\partial _b \ldots \partial _c \phi ^a$
and
$\partial ^k \phi ^a = \partial _m\partial _p \ldots \partial _q \phi ^a$
products of
$l$
and
$k$
derivatives, respectively. The terms involving
$\phi ^a$
are not part of the integral but the fact that
$\phi ^a$
is harmonic plays a key role in finding a simple result.
The first point to note is that, if the total number of normals
$l+k+1$
is odd, the integral vanishes because it changes sign upon the reflection
$\boldsymbol {x}\rightarrow -\boldsymbol {x}$
. Hence,
$l+k+1=2N$
, with
$N$
an integer. It is known that, over the unit sphere, such an integral equals
$4\pi /(2N+1)!!$
multiplied by the sum of
$(2N-1)!!$
products of
$N$
Kronecker
$\delta$
values, each one carrying a pair of unequal indices chosen from the set
$(a,\,b,\ldots ,c),\, (m,\, p,\ldots q)$
and
$i$
(Google AI Mode 2026). It is evident that any such product containing even a single
$\delta$
with both indices from the same set, either
$(a,\,b,\ldots ,c)$
or
$(m,\, p,\ldots q)$
, such as
$\delta _{ab}$
or
$\delta _{pq}$
, will vanish as it will include the term
${\nabla} ^2\phi ^a$
. Thus the only non-zero contributions can come from products of
$\delta$
values in which one index, chosen from the
$l+1$
indices
$(i,\,a,\,b,\ldots ,c)$
, is paired with one index chosen from the
$k$
indices
$(m,\, p,\ldots q)$
. For such non-zero contributions to exist the number of indices in the two sets must be the same as, otherwise, two additional indices in one set must necessarily introduce a
$\delta$
with both indices from that set. Since the index
$i$
is prescribed, the only non-zero contributions must have either
$k=l-1$
or
$k=l+1$
. On the basis of these considerations (3.5) becomes
\begin{align} \frac {1}{2}\oint _S (u^a)^2_i {\rm d}S_x = \frac {1}{2} \left (d_0d_1 I_i^{(0,1)} + d_1 (d_0 I_i^{(1,0)}+d_2 I_i^{(1,2)} )+ d_2 d_1I_i^{(2,1)} + d_3I_i^{(2,3}) +\ldots \right ) , \end{align}
which, due to the symmetry of
$I_i^{(l,k)}$
, simply becomes
$d_0 d_1I_i^{(0,1)} +d_1d_2 I_i^{(1,2)}+\ldots$
, as in (3.6).
It is sufficient, therefore, to limit our considerations to integrals of the form
$ I_i^{(l-1,l)}$
which, on the basis of the previous considerations, once contracted with the derivatives of
$\phi ^a$
as in (A1), will have the same effect as
\begin{align} \frac {4\pi R^2}{(2l+1)!!}C_l \delta _{im}\delta _{ap} \ldots \delta _{cq} , \end{align}
with suitable constants
$C_l$
. To determine these constants we set
$n_a=r_a/R$
etc. and re-write
$ I^{(l-1,l)}_i$
as
\begin{align} I^{(l-1,l)}_i = R^{-(2l-1)}\oint _S \underbrace {\left (r_ar_b \ldots r_c\right )}_{l-1 \,\,{ terms}}\underbrace {\left (r_m r_p \ldots r_q\right )}_{l\,\,{terms}}n_i {\rm d}S_x \, \big(\partial ^{l-1}\phi ^a \big) \big(\partial ^l \phi ^a \big) . \end{align}
Upon using the divergence theorem to convert the surface integral to an integral over the volume
$V$
of the particle, by the symmetry with respect to the sets of indices
$(a,\,b,\ldots ,\,c)$
and
$(m,\,p,\ldots ,\,q)$
of the integrand, we may write the result as
\begin{eqnarray} I^{(l-1,l)}_i &=& (l-1)R^{-(2l-1)}\delta _{ia}\int _V \underbrace {\left (r_b \ldots r_c\right )}_{l-2 \,\,{terms}}\underbrace {\left (r_m r_p \ldots r_q\right )}_{l\,\,{terms}}{\rm d}V \, \big(\partial ^{l-1}\phi ^a \big) \big(\partial ^l \phi ^a \big) \nonumber \\&&\qquad + lR^{-(2l-1)}\delta _{ip}\int _V \underbrace {\left (r_ar_b \ldots r_c\right )}_{l-1 \,\,{terms}}\underbrace {\left (r_p \ldots r_q\right )}_{l-1\,\,{terms}}{\rm d}V \, \big(\partial ^{l-1}\phi ^a \big) \big(\partial ^l \phi ^a \big) \, . \end{eqnarray}
The first term must vanish because each one of its
$\delta$
-products must also contain a
$\delta$
with both indices from the second group. For the second term we note that the volume integral over a sphere of radius
$R$
can be written as
$ \oint (\ldots ){\rm d}V = \int _0^R r^2 {\rm d}r \oint (\ldots ) {\rm d}\varOmega$
, with
$\varOmega$
the solid angle. Thus
\begin{eqnarray} I^{(l-1,l)}_i &=& lR^{-(2l-1)}\delta _{im}\int _0^R r^2 r^{2l-2} {\rm d}r\oint \underbrace {\left (n_an_b \ldots n_c\right )}_{l-1 \,\,{terms}}\underbrace {\left (n_p \ldots n_q\right )}_{l-1\,\,{terms}} {\rm d}\varOmega \big(\partial ^{l-1}\phi ^a \big) \big(\partial ^l\phi ^a \big) \nonumber \\&=& \frac {l}{2l+1}R^2\delta _{im}\oint \underbrace {\left (n_b \ldots n_c\right )}_{l-2 \,\,{ terms}}\underbrace {\left (n_p \ldots n_q\right )}_{l-1\,\,{terms}} n_a {\rm d}\varOmega \big(\partial ^{l-1}\phi ^a \big) \big(\partial ^l\phi ^a \big) . \end{eqnarray}
Using (A3), with the sphere’s radius taken as
$R$
on the left-hand side and as 1 on the right-hand side, we find
\begin{eqnarray} \frac {4\pi R^2}{(2l+1)!!}C_l \delta _{im}\delta _{ap}\ldots \delta _{cq} = \frac {l}{2l+1}R^2 \delta _{im}\left [ \frac {4\pi }{(2l-1)!!}C_{l-1}\delta _{ap} \ldots \delta _{cq}\right ] \, . \end{eqnarray}
Hence, we deduce that
$C_l=lC_{l-1}$
so that
$C_l=l!C_1$
. The integral
$I^{(0,1)}_i$
is elementary and, upon comparison with (A3), shows that
$C_1=1$
so that
$C_l=l!$
. We therefore conclude that
\begin{align} I^{(l-1,l)}_i = \frac {4 l!}{(2l+1)!!}\pi R^2 (\underbrace {\partial _a\partial _b \ldots \partial _c}_{l-1 \,\, {terms}} \phi ^a)(\partial _i \underbrace {\partial _a\partial _b \ldots \partial _c}_{l-1 \,\, {terms}}\phi ^a) , \end{align}
from which both (3.7) and (4.6) follow. A slightly different form of the integral arises in the process of calculating the kinetic energy. It is easy to see that this form is readily reducible to that for
$I_i^{(l-1,l)}$
.
