On Spheres and Spheroids in Stokes Flow - 8/26
Introduction
I spent a good amount of this summer combing through Kim and Karilla’s Microhydrodynamics. Turns out you can do really very cool things when the governing equations of your problem (Stokes equations) are linear. A few months ago, I had no idea the rabbit hole was so deep.
I will almost certainly be doing the bulk of my PhD work on things related to this, so I figured I should write up some half-decent notes. And by half-decent, I mean the interesting bits of the first dozen chapters of Microhydrodynamics.
This is a walkthrough of the singularity system for spheres, then for the much more complex ellipsoids, then for prolate spheroids, and then for slender bodies. There are some applications throughout and discussed more at the end. These notes are a bit rougher than other ones, so proceed with caution.
Singularity Systems
The Green’s function is our most fundamental tool for studying solutions to the Stokes Equations. Once we develop the singularity solution and integral representations, we may proceed to spherical and spheroid solutions.
The Stokes Equations
We begin with the Stokes equations, \begin{equation} -\partial_i p\left(\mathbf{x}\right) +\mu\partial^2 u_i\left(\mathbf{x}\right) = -f_i\left(\mathbf{x}\right),\quad\partial_iu_i\left(\mathbf{x}\right)=0. \end{equation} \(u_i\left(\mathbf{x}\right)\) is the velocity field, \(p\left(\mathbf{x}\right)\) is the pressure field, and \(f_i\) is a non-conservative force distribution. \(\mu\) is the fluid viscosity. The equation can also be written as the divergence of a stress tensor, \begin{equation} \partial_j\sigma_{ij}\left(\mathbf{x}\right) = -f_i\left(\mathbf{x}\right),\quad\sigma_{ij}\left(\mathbf{x}\right)=-p\left(\mathbf{x}\right)\delta_{ij}+2\mu e_{ij}\left(\mathbf{x}\right), \end{equation} where \(e_{ij}\) are the components of the symmetric rate of deformation tensor, \begin{equation} e_{ij}\left(\mathbf{x}\right)=\frac{1}{2}\left(\partial_i u_j\left(\mathbf{x}\right) + \partial_j u_i\left(\mathbf{x}\right)\right). \end{equation}
The Green’s Function for Stokes Flow
We wish to solve the Stokes Equations in the presence of a Dirac delta function forcing, \begin{equation}\label{eq: delta function stokes} -\partial_i p\left(\mathbf{x}\right) +\mu\partial^2 u_i\left(\mathbf{x}\right) = -F_i\delta\left(\mathbf{x}\right)\quad\partial_iu_i=0. \end{equation} to obtain a Green’s function solution for our flow. This is most straightforward to do in wavenumber space. Taking a Fourier transform of Eq. \ref{eq: delta function stokes}, the equation is evaluated in \(k\)-space, \begin{equation} -i k_i \tilde{p} -\mu k^2 \tilde{u}_i = -F_i\quad k_i\tilde{u}_i=0. \end{equation} It is not difficult to obtain the solution of the linear system, \begin{equation} \tilde{u}_i = \frac{F_j}{\mu k^2}\left(\delta _{ij}-\frac{k _ik _j}{k^2}\right), \quad \tilde{p} = -\frac{i k_i F_i}{k^2}. \end{equation} Reverse Fourier-transforming, we return to position space with, \begin{equation} u_i = \frac{F _j}{8\pi\mu}\mathcal{G} _{ij}, \quad p = \frac{F_i}{8\pi \mu}\mathcal{P} _i, \end{equation} where \(\mathcal{G} _{ij}\) is the Oseen Tensor and \(\mathcal{P}_i\) is the corresponding pressure field tensor, given by, \begin{equation} \mathcal{G} _{ij} = \frac{1}{r}\delta _{ij}+\frac{1}{r^3}x_ix_j, \end{equation} \begin{equation} \mathcal{P} _i = 2\mu\frac{x _i}{r^3}. \end{equation} Now we may also construct the stress field tensor, \begin{equation} \sigma _{ij}=\frac{F_k}{8\pi \mu}\left(-\mathcal{P} _k\delta _{ij} + \mu\left(\partial_i\mathcal{G} _{jk}+\partial_j\mathcal{G} _{ik}\right)\right) = F_k\Sigma _{ijk}, \end{equation} with, \begin{equation} 8\pi\mu\Sigma _{ijk} = -\mathcal{P} _k\delta _{ij} + \mu\left(\partial_i\mathcal{G} _{jk}+\partial_j\mathcal{G} _{ik}\right) = -\mu\frac{6}{r^5}x_ix_jx_k. \end{equation} We can actually remove the force term to rewrite Eq. \ref{eq: delta function stokes} in the more general form, placing a singularity at \(\mathbf{\xi}\), \begin{equation}\label{eq: delta function stokes tensor representation 1} -\partial_i \mathcal{P} _k\left(\mathbf{x}-\mathbf{\xi}\right) +\mu\partial^2 \mathcal{G} _{ik}\left(\mathbf{x}-\mathbf{\xi}\right) = 8\pi\mu\partial _j\Sigma _{ijk}\left(\mathbf{x} - \mathbf{\xi}\right) = -8\pi\mu \delta _{ik} \delta\left(\mathbf{x}-\mathbf{\xi}\right). \end{equation} \begin{equation}\label{eq: delta function stokes tensor representation 2} \partial_i\mathcal{G} _{ij}\left(\mathbf{x} - \mathbf{\xi}\right)=0. \end{equation}
The Integral Representation with Green’s Functions
We derived the Green’s functions to power a representation of the velocity field as a surface integral. First, we will need the reciprocal theorem.
Consider a closed region of fluid \(V\) enclosed by a surface \(S\). Two velocity fields \(u_i^1\) and \(u_i^2\) satisfy the Stokes Equations with different boundary conditions on the surface \(S\). Now notice, \begin{equation} \sigma^1 _{ij} e^2 _{ij} = 2\mu e^1 _{ij}e^2 _{ij} = \sigma^2 _{ij} e^1 _{ij}. \end{equation} But consider, \begin{equation} \sigma^1 _{ij}e^2 _{ij}=\sigma^1 _{ij}\partial _j u^2_i = \partial _j\left(\sigma^1 _{ij}u^2_i\right)-u^2 _i\partial _j\sigma^1 _{ij}=\partial _j\left(\sigma^2 _{ij}u^1_i\right)-u^1 _i\partial _j\sigma^2 _{ij}. \end{equation} And integrating over the region \(V\) and applying the divergence theorem, where \(n_i\) is a unit vecotr pointing out of the fluid-filled region \(V\), \begin{equation}\label{eq: Reciprocal theorem} \oint _S dS\ \sigma^1 _{ij}u^2_i n_j - \int_V dV\ u^2 _i\partial _j\sigma^1 _{ij} = \oint _S dS\ \sigma^2 _{ij}u^1_i n_j - \int_V dV\ u^1 _i\partial _j\sigma^2 _{ij}. \end{equation} Eq. \ref{eq: Reciprocal theorem} is known as the Reciprocal theorem.
We can also apply the Reciprocal theorem to the Green’s function representation. Let \(u^1_i\) correspond to the solution of the Stokes Equations in the presence of a point forcing at \(\mathbf{x}\) and let \(u^2\) be replaced by \(u\), the velocity field in the region of interest. Eq. \ref{eq: Reciprocal theorem} becomes, \begin{equation} \oint _S dS\left(\mathbf{\xi}\right)\ F_k\Sigma _{ijk}\left(\mathbf{\xi}-\mathbf{x}\right)u_i\left(\mathbf{\xi}\right) n_j\left(\mathbf{\xi}\right)+\int_V dV\left(\mathbf{\xi}\right)\ u _i\left(\mathbf{\xi}\right) F_i\delta\left(\mathbf{\xi}-\mathbf{x}\right)\nonumber \end{equation} \begin{equation} = \oint _S dS\left(\mathbf{\xi}\right)\ \sigma _{ij}\left(\mathbf{\xi}\right)\frac{F_k}{8\pi \mu}\mathcal{G} _{ik}\left(\mathbf{\xi}-\mathbf{x}\right)n_j\left(\mathbf{\xi}\right) - \int_V dV\left(\mathbf{\xi}\right)\ \frac{F_k}{8\pi \mu}\mathcal{G} _{ik}\left(\mathbf{\xi}-\mathbf{x}\right)\partial _j\sigma _{ij}\left(\mathbf{\xi}\right). \end{equation}
As this holds for any \(F_k\) forcing, we can do away with this index, and assuming the \(u\) field is force-free such that \(\partial _j\sigma _{ij}=0\), \begin{equation} \frac{1}{8\pi\mu}\oint _S dS\left(\mathbf{\xi}\right)\ \sigma _{ij}\left(\mathbf{\xi}\right)\mathcal{G} _{ik}\left(\mathbf{\xi}-\mathbf{x}\right)n_j\left(\mathbf{\xi}\right) =\oint _S dS\left(\mathbf{\xi}\right)\ \Sigma _{ijk}\left(\mathbf{\xi}-\mathbf{x}\right)u_i\left(\mathbf{\xi}\right) n_j\left(\mathbf{\xi}\right)+\int_V dV\left(\mathbf{\xi}\right)\ u _k\left(\mathbf{\xi}\right) \delta\left(\mathbf{\xi}-\mathbf{x}\right). \end{equation} Now we can exchange the order of \(\mathbf{x}\) and \(\mathbf{\xi}\) using \(\mathcal{G} _{ij}\left(-\mathbf{x}\right)=\mathcal{G} _{ij}\left(\mathbf{x}\right)\) and \(\Sigma _{ijk}\left(-\mathbf{x}\right)=-\Sigma _{ijk}\left(\mathbf{x}\right)\). For the velocity field inside the region \(V\), the \(\delta\) function integral turns into \(u_k\left(\mathbf{x}\right)\), while outside this region, it evaluate to \(0\), giving the following surface integral representation for the velocity, \begin{equation}\label{eq: surface integral velocity representation 1} \text{if}\ \mathbf{x}\in V,\quad u_k\left(\mathbf{x}\right)=\frac{1}{8\pi\mu}\oint _S dS\left(\mathbf{\xi}\right)\ \sigma _{ij}\left(\mathbf{\xi}\right)\mathcal{G} _{ik}\left(\mathbf{x}-\mathbf{\xi}\right)n_j\left(\mathbf{\xi}\right) + \oint _S dS\left(\mathbf{\xi}\right)\ \Sigma _{ijk}\left(\mathbf{\xi}-\mathbf{x}\right)u_i\left(\mathbf{\xi}\right) n_j\left(\mathbf{\xi}\right), \end{equation} \begin{equation}\label{eq: surface integral velocity representation 2} \text{if}\ \mathbf{x}\notin V,\quad 0=\frac{1}{8\pi\mu}\oint _S dS\left(\mathbf{\xi}\right)\ \sigma _{ij}\left(\mathbf{\xi}\right)\mathcal{G} _{ik}\left(\mathbf{x}-\mathbf{\xi}\right)n_j\left(\mathbf{\xi}\right) + \oint _S dS\left(\mathbf{\xi}\right)\ \Sigma _{ijk}\left(\mathbf{\xi}-\mathbf{x}\right)u_i\left(\mathbf{\xi}\right) n_j\left(\mathbf{\xi}\right). \end{equation}
This is helpful, but we are also interested in the case of a particle immersed in a fluid region. Consider now a fluid region \(V_f\) bounded by a particle surface \(S_p\) and a far-away surface \(S_\infty\). The particle undergoes some rigid-body motion such that \(u_i=u^p_i\) on the surface of and inside the particle, and far away from the particle there is ambient flow \(u_i=u^\infty_i\) on \(S_\infty\). We are interested in evaluating the disturbance field \(u^D_i=u_i-u^\infty_i\) in the region \(V_f\). With Eq. \ref{eq: surface integral velocity representation 1} and Eq. \ref{eq: surface integral velocity representation 2}, we see, \begin{equation} \text{if}\ \mathbf{x}\in V_f,\quad u^D_k\left(\mathbf{x}\right)=-\frac{1}{8\pi\mu}\oint _{S_p} dS\left(\mathbf{\xi}\right)\ \sigma^D _{ij}\left(\mathbf{\xi}\right)\mathcal{G} _{ik}\left(\mathbf{x}-\mathbf{\xi}\right)\hat{n}_j\left(\mathbf{\xi}\right)-\oint _{S _p} dS\left(\mathbf{\xi}\right)\ \left(u^p_i-u^\infty_i\right)\Sigma _{ijk}\left(\mathbf{\xi}-\mathbf{x}\right) \hat{n}_j\left(\mathbf{\xi}\right)\nonumber \end{equation} \begin{equation}\label{eq: integral represntation outside particle} -\frac{1}{8\pi\mu}\oint _{S _\infty} dS\left(\mathbf{\xi}\right)\ \sigma^D _{ij}\left(\mathbf{\xi}\right)\mathcal{G} _{ik}\left(\mathbf{x}-\mathbf{\xi}\right)\hat{n}_j\left(\mathbf{\xi}\right) -\oint _{S _\infty} dS\left(\mathbf{\xi}\right)\ u^D_i\Sigma _{ijk}\left(\mathbf{\xi}-\mathbf{x}\right) \hat{n}_j\left(\mathbf{\xi}\right), \end{equation} where \(\hat{n}_i\) points into the fluid region in the direction opposite to \(n_i\). The integrals over \(S _\infty\) vanish as the edge of the boundary approaches infinity if \(u^D_i\rightarrow0\) and \(r\sigma^D _{ij}\rightarrow0\) as \(r\rightarrow\infty\) as \(\mathcal{G} _{ij}\sim r^{-1}\) and \(\Sigma _{ijk}\sim r^{-2}\) and \(dS\sim r^2\). The first condition holds, \(u^D_i\rightarrow0\), as we expect the influence of the particle to vanish far from its surface and by the definition of \(u_i^\infty\). We similarly expect \(\sigma^D _{ij}\sim r^{-2}\) since the gradient of the disturbance field goes with, \(\partial_i u^D_j\rightarrow r^{-2}\).
Now apply the same integral representation to the region inside the particle where the velocity \begin{equation} \text{if}\ \mathbf{x}\in V _p, \quad u_k^p- u^\infty_k = \frac{1}{8\pi\mu}\oint _{S_p} dS\left(\mathbf{\xi}\right)\ \left(\sigma^p _{ij}-\sigma^\infty _{ij}\left(\mathbf{\xi}\right)\right)\mathcal{G} _{ik}\left(\mathbf{x}-\mathbf{\xi}\right)\hat{n}_j\left(\mathbf{\xi}\right)\nonumber \end{equation} \begin{equation}\label{eq: integral represntation inside particle} +\oint _{S _p} dS\left(\mathbf{\xi}\right)\ \left(u^p_i-u^\infty_i\right)\Sigma _{ijk}\left(\mathbf{\xi}-\mathbf{x}\right) \hat{n}_j\left(\mathbf{\xi}\right). \end{equation}
The surface traction integral, \(\sigma^p _{ij}\) vanishes, since its net contribution is a constant pressure term in the case of a rigid body and the integral over the solenoidal Green’s function vanishes.
We add Eq. \ref{eq: integral represntation outside particle} and Eq. \ref{eq: integral represntation inside particle} to obtain a representation of the velocity field everywhere. Integrals over \(S_\infty\) are also removed since they vanish as discussed above, \begin{equation}\label{eq: greens function integral rep} u_k\left(\mathbf{x}\right) = u_k^\infty -\frac{1}{8\pi\mu}\oint _{S_p} dS\left(\mathbf{\xi}\right)\ \sigma _{ij}\left(\mathbf{\xi}\right)\mathcal{G} _{ik}\left(\mathbf{x}-\mathbf{\xi}\right)\hat{n}_j\left(\mathbf{\xi}\right). \end{equation} So the disturbance field everywhere in the fluid domain can be expressed as a surface integral with a Green’s function kernel over the stress at the particle surface.
The Multipole Expansion
Now we may expand the integral representation in the far field, the limit of \(|\mathbf{x}|\gg|\mathbf{\xi}|\) so that \(\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{\xi}\right)\approx\mathcal{G} _{ij}\left(\mathbf{x}\right)\). Formally, expanding the Green’s function around \(\mathbf{\xi}=0\),
\begin{equation} \mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{\xi}\right)=\sum _{n=0}^\infty\frac{ \left(-1\right)^n }{n!}\left(\xi _k\partial _k\right)^n\mathcal{G} _{ij}\left(\mathbf{x}\right). \end{equation}
Inserting this into the velocity integral,
\begin{equation} u_k\left(\mathbf{x}\right) -u_k^\infty = -\frac{1}{8\pi\mu}\sum _{n=0}^\infty\frac{ \left(-1\right)^n }{n!}\oint _{S_p} dS\left(\mathbf{\xi}\right)\ \left(\sigma _{ij} \hat{n}_j\left(\mathbf{\xi}\right) \xi _{k_1}\xi _{k_2}…\xi _{k_n}\right)\partial _{k_1}\partial _{k_2}…\partial _{k_n}\mathcal{G} _{ik}\left(\mathbf{x}\right), \end{equation}
\begin{equation} u _k\left(\mathbf{x} \right) -u _k^\infty = -F _i\frac{\mathcal{G} _{ik}\left(\mathbf{x}\right)}{8\pi\mu} +D _{ij}\frac{\partial_j\mathcal{G} _{ik}\left(\mathbf{x}\right)}{8\pi\mu} + …, \end{equation}
with, \begin{equation} F _i = \oint _{S_p} dS \left(\mathbf{\xi}\right)\ \sigma _{ij} \hat{n}_j \left( \mathbf{\xi} \right), \end{equation}
\begin{equation} D _{ij} = \oint _{S_p} dS\left(\mathbf{\xi}\right)\ \sigma _{il} \hat{n}_l\left(\mathbf{\xi}\right)\xi _j . \end{equation}
\(F _i\) is known as the stokeslet, the next effect is the dipole \(D _{ij}\) that decays with \(r^{-2}\). Because the isotropic portion of \(D _{ij}\) has no dynamic significance since \(\partial_i\mathcal{G} _{ij}=0\), we can split the term into a symmetric stresslet and an antisymmetric rotlet, \begin{equation} D _{ij} -\frac{1}{3}D _{kk}\delta _{ij} = S _{ij} + T _{ij}, \end{equation} \begin{equation} S _{ij} = \frac{1}{2}\oint _{S_p} dS\left(\mathbf{\xi}\right)\ \left(\sigma _{il}\xi _j + \sigma _{jl} \xi _i\right)\hat{n}_l\left(\mathbf{\xi}\right) - \frac{1}{3}\oint _{S_p} dS\left(\mathbf{\xi}\right)\ \sigma _{lm}\xi _l \hat{n}_m\left(\mathbf{\xi}\right)\delta _{ij}, \end{equation} \begin{equation} T _{ij} = \frac{1}{2}\oint _{S_p} dS\left(\mathbf{\xi}\right)\ \left(\sigma _{il}\xi _j - \sigma _{jl} \xi _i\right)\hat{n}_l\left(\mathbf{\xi}\right). \end{equation} We may also identify the relation of the rotlet with the torque pseudovector, \begin{equation} -\epsilon _{ijk}T _{jk} = -\epsilon _{ijk}\oint _{S_p} dS\left(\mathbf{\xi}\right)\ \sigma _{jl} \xi _k \hat{n}_l\left(\mathbf{\xi}\right) = T _i, \end{equation} and, \begin{equation} T _{jk}=-\frac{1}{2}\epsilon _{jkl}T _l. \end{equation} Next, we will use the multipole expansion to identify the velocity field around a spherical particle
The Singularity System for Spheres
First, we calculate a few derivatives of the Oseen tensor, which will come in handy, \begin{equation} \mathcal{G} _{ij} = \frac{1}{r}\delta _{ij}+\frac{1}{r^3}x_ix_j, \end{equation} \begin{equation} \partial_k\mathcal{G} _{ij} = \frac{1}{r^3}\left(-\delta _{ij}x_k+\delta _{ik}x_j+\delta _{jk}x_i\right)-\frac{3}{r^5}x _ix _jx _k, \end{equation} \begin{equation} \partial^2\mathcal{G} _{ij} = \frac{2}{r^3}\delta _{ij}-\frac{6}{r^5}x_ix_j, \end{equation} \begin{equation} \partial_k\partial^2\mathcal{G} _{ij} = -\frac{6}{r^5}\left(\delta _{ij}x_k+\delta _{ik}x_j+\delta _{jk}x_i\right)+\frac{30}{r^7}x _ix _jx _k. \end{equation} \begin{equation} \partial^2\partial^2\mathcal{G} _{ij} = 0. \end{equation} This last identity is expected since the velocity field is biharmonic by definition.
We consider the velocity prescribed on the spherical surface located at the origin as a superposition of translational, rotational, and symmetric \(\left(E_{ij}=E_{ji}\right)\) straining modes, \begin{equation} u^p_i = U_i + \epsilon_{ijk}\Omega_jx_k + E_{ij}x_j. \end{equation} Generally, we may also explore higher-order components, but these are calculated in much the same way and will end up corresponding to higher-order modes of the multipole expansion, which are of less interest to us. The sphere also sits in an ambient flow, with the far-field velocity taking the form, \begin{equation} u^\infty_i = U^\infty_i + \epsilon_{ijk}\Omega^\infty_jx_k + E^\infty_{ij}x_j, \end{equation} where \(U^\infty_i\), \(\Omega^\infty_j\), and \(E^\infty_{ij}\) are constant. The boundary conditions describe the velocity at the sphere’s surface, \begin{equation} u_i\left(r=a\right)=u^p_i \end{equation} where \(a\) is the radius of the sphere, and the velocity far from the sphere takes the far-field unperturbed values, \begin{equation} u_i\left(r\rightarrow\infty\right)=u^\infty_i. \end{equation} Note that we have assumed we know the velocity of the sphere within the moving fluid. Suppose instead we know the velocity of the sphere in a quiescent fluid and are interested in the dynamics in the presence of \(u^\infty_i\). This could be the case if we are studying, for instance, a swimming body. In this case, if the body would translate with an intrinsic velocity \(u^s_i\) in the absence of \(u^\infty_i\), we simply replace \(u^p_i\rightarrow u^s_i+u^\infty_i\).
Considering first the translation case. Notice at \(r=a\), \begin{equation} \left(1+\frac{a^2}{6}\partial^2\right)\mathcal{G} _{ij}\left(r=a\right) = \frac{4}{3a}\delta _{ij}, \end{equation} from which the velocity is represented as, \begin{equation} u_i\left(\mathbf{x}\right)=-6\pi\mu a\left(U^\infty_j-U_j\right)\left(1+\frac{a^2}{6}\partial^2\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu} + U^\infty_i, \end{equation} for a translating sphere, the multipole expansion has only these two terms, and only two derivative terms of the Green’s function are necessary.
Now for the rotating case. At \(r=a\), \begin{equation} \epsilon_{jkl}\epsilon_{nml}\partial_m\mathcal{G} _{in}\left(r=a\right) = \partial_k\mathcal{G} _{ij}\left(r=a\right)-\partial_j\mathcal{G} _{ik}\left(r=a\right) = \frac{2}{a^3}\left(\delta _{ik}x_j-\delta _{ij}x_k\right) = \frac{2}{a^3}\epsilon _{jkl}\epsilon _{ilm}x_m. \end{equation} So set the velocity field, \begin{equation} u_i\left(\mathbf{x}\right)=4\mu\pi a^3\epsilon _{ljk}\frac{\partial_l\mathcal{G} _{ik}\left(\mathbf{x}\right)}{8\pi\mu}\left(\Omega_j-\Omega^\infty_j\right) + \epsilon _{ijk}\Omega^\infty _jx _k \end{equation}
And lastly, for the strained case. At \(r=a\), \begin{equation} \left(1+\frac{a^2}{10}\partial^2\right)\left(\partial_k\mathcal{G} _{ij}\left(r=a\right)+\partial_j\mathcal{G} _{ik}\left(r=a\right)\right) = -\frac{6}{5a^3}\left(\delta _{ij}x_k+\delta _{ik}x_j-\frac{2}{3}\delta _{jk}x_i\right), \end{equation} As \(E _{ij}=E _{ji}\) and \(E _{kk}=0\), \begin{equation} E _{jk}\left(1+\frac{a^2}{10}\partial^2\right)\partial_k\mathcal{G} _{ij}\left(r=a\right) = \frac{1}{2}E _{jk}\left(1+\frac{a^2}{10}\partial^2\right)\left(\partial_k\mathcal{G} _{ij}\left(r=a\right)+\partial_j\mathcal{G} _{ik}\left(r=a\right)\right) = -\frac{6}{5a^3}E _{ik}x_k, \end{equation} and we identify the velocity field as, \begin{equation} u_i\left(\mathbf{x}\right)= \frac{20}{3}\pi\mu a^3 \left(E _{jk}^\infty-E _{jk}\right) \left(1+\frac{a^2}{10}\partial^2\right)\frac{\partial_k\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu} + E^\infty _{ij}x _j. \end{equation}
In summary, we may express our velocity field as a sum of the translational, rotational, and straining modes, \begin{equation} u_i\left(\mathbf{x}\right)= \left[-F_j\left(1+\frac{a^2}{6}\partial^2\right) + S_{jk}\left(1+\frac{a^2}{10}\partial^2\right)\partial_k -\frac{1}{2} \epsilon_{jkl}T _k\partial_l\right]\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu} + u^\infty _i, \end{equation} \begin{equation}\label{eq: force translating particle const background} F_i = 6\pi\mu a\left(U^\infty _i -U _i\right), \end{equation}
\begin{equation} S_{jk} = \frac{20}{3}\pi\mu a^3 \left(E _{jk}^\infty-E _{jk}\right), \end{equation} \begin{equation} T _k = 8\mu\pi a^3\left( \Omega^\infty_k-\Omega_k\right). \end{equation} \(F _i\), \(S _{jk}\), \(T _k\) are the hydrodynamic force, stress, and torque associated with the particle’s motion, as in this is the force exerted by the fluid onto the particle. Simply reverse the signs for the force experienced by the fluid.
Faxén Laws
We require one last element, that being the inverse relations. We know how to find the flow given the background fluid motion and the behaviour of the sphere, but suppose we are instead given a non-constant background flow \(u^\infty_i\left(\mathbf{x}\right)\), how do we determine the moments \(F_i\left(\mathbf{x}\right)\), \(T_i\left(\mathbf{x}\right)\), and \(S_{ij}\left(\mathbf{x}\right)\)?
We first consider purely translational motion, starting again with the reciprocal theorem, Eq. \ref{eq: Reciprocal theorem}, \begin{equation}\label{eq: flaxen reciprocal theorem} \oint _{S_p} dS\ \sigma^1 _{ij}u^2_i \hat{n}_j + \int_V dV\ u^2 _i\partial _j\sigma^1 _{ij} = \oint _{S_p} dS\ \sigma^2 _{ij}u^1 _i \hat{n} _j + \int_V dV\ u^1 _i\partial _j\sigma^2 _{ij}, \end{equation} where \(\hat{n}_i\) is the vector pointing out of the particle into the fluid. Now let \(u^1 _i\) be the solution for a spherical particle translating with velocity \(U_i\) in a quiescent fluid. Note that this means that \(\partial _j\sigma _{ij}^1=0\) since the auxiliary flow satisfies the homogeneous Stokes momentum equation in the fluid. For the second velocity field, \(u^2 _i\), choose the velocity field generated by a point force \(P _i\delta\left(\mathbf{x}-\mathbf{y}\right)\) where \(\mathbf{y}\) lies outside of the particle and the particle is stationary. Because the second particle is stationary, we have \(u_i^2=0\) on the particle surface, and \(\partial _j\sigma _{ij}^2=-P _i\delta\left(\mathbf{x}-\mathbf{y}\right)\) from the point force. Inserting these fields into Eq. \ref{eq: flaxen reciprocal theorem},
\begin{equation} F^2_i U_i - P_iu^1 _i\left(\mathbf{y}\right) = 0, \end{equation}
where \(F^2_i\) is the hydrodynamic force on the second particle,
\begin{equation} F^2_i = \oint _{S_p} dS\ \sigma^2 _{ij} \hat{n} _j. \end{equation}
We can write \(u^1_i\) with our singularity solution,
\begin{equation} u^1 _i = U _j \mathcal{F}\left[ \frac{\mathcal{G} _{ij}\left( \mathbf{x}-\mathbf{\xi} \right) }{8\pi\mu} \right], \end{equation} \begin{equation} \mathcal{F}\left[f\right]=6\pi\mu a\left( 1+\frac{a^2 }{6}\partial^2 \right)f, \end{equation}
where \(\mathcal{F}\) is a linear functional and \(\mathbf{\xi}\) is the centre of the sphere, from which we identify
\begin{equation} F^2_i = P_j \mathcal{F} \left[ \frac{\mathcal{G} _{ij}\left(\mathbf{y}-\mathbf{\xi}\right)}{8\pi\mu} \right]. \end{equation}
But the element in the functional is nothing more than the ambient \(u_i^2\) field evaluated at the sphere centre,
\begin{equation} P_j \mathcal{F} \left[ \frac{\mathcal{G}_{ij} \left(\mathbf{y}-\mathbf{\xi}\right) }{8\pi\mu} \right] = \mathcal{F} \left[u^2_i\left( \mathbf{x} \right)\right] _{\mathbf{x}=\mathbf{\xi}}. \end{equation}
Consider now an ambient flow field \(u_i^\infty\left(\mathbf{x}\right)\). This flow field can be generated using our Green’s function integral representation, Eq. \ref{eq: greens function integral rep}. So if the second particle sits in an ambient field \(u_i^\infty\left(\mathbf{x}\right)\) constructed from a superposition of Green’s functions, we may identify the force on the second particle as,
\begin{equation}\label{eq: force flaxen relation} F_i = \mathcal{F} \left[ u_i^\infty\left(\mathbf{x}\right) \right]_{\mathbf{x}=\mathbf{\xi}} = 6\pi\mu a\left( 1+\frac{a^2 }{6}\partial^2 \right)u_i^\infty\left(\mathbf{x}\right)| _{\mathbf{\xi}}. \end{equation}
The force on the stationary particle is obtained using the same functional that gave us our multipole expansions. Notice that if the ambient flow is constant, \(u_i^\infty\left(\mathbf{x}\right)=U^\infty_i\), we recover Eq. \ref{eq: force translating particle const background}.
By a similar argument with the reciprocal theorem, one can find the rotlet and stresslet, \begin{equation} T_i = 4\pi\mu a^3 \epsilon_{ijk}\partial_ju_k^\infty| _{\mathbf{\xi}}, \end{equation}
\begin{equation}\label{eq: stresslet flaxen relation} S_{ij} = \frac{20}{3}\pi\mu a^3\left( 1+\frac{a^2 }{10}\partial^2 \right)\frac{1}{2}\left(\partial_iu_j^\infty+\partial_ju_i^\infty\right)| _{\mathbf{\xi}}. \end{equation}
Example Motion for Rigid Body Swimming
This is somewhat supplementary, but helpful in understanding our Faxén and multipole relations. Suppose we have a spherical swimmer in a flow \(u_i^\infty\left(\mathbf{x}\right)\). In the absence of the flow, the swimmer would undergo rigid body motion with velocity \(U^s_i\) and angular velocity \(\Omega^s_i\), so that the velocity of the flow at the sphere’s surface is given by, \begin{equation} u_i\left(r=a\right)=U^s_i+\epsilon_{ijk}\Omega^s_jr_k+u_i^\infty\left(r=a\right). \end{equation} We are interested in the resulting translational and rotational velocity of our swimmer in the flow, denoted by \(U^p_i\) and \(\Omega^p_i\). From our Faxén relation, Eq. \ref{eq: force flaxen relation}, the hydrodynamic force on the sphere is, \begin{equation} F^\text{hydro}_i = 6\pi\mu a\left( 1+\frac{a^2 }{6}\partial^2 \right)\left(u_i^\infty\left(\mathbf{x}\right)-U^p_i\right)| _{\mathbf{\xi}}. \end{equation} The sphere also experiences a propulsive swimming force \(F_i^\text{swim}\), producing the translational swimming with velocity \(U^s_i\). This force is found with Eq. \ref{eq: force translating particle const background}, \begin{equation} F^\text{swim}_i = 6\pi\mu a U^s_i, \end{equation} where the sign has changed since \(F_i^\text{swim}\) is the force required by the sphere to move through a quiescent fluid. The sphere’s internal propulsive force must balance with the next hydrodynamic force, as there is no net force experienced by the swimming particle, resulting in, \begin{equation} F^\text{swim}_i + F^\text{hydro}_i = 0. \end{equation} From this, it is not difficult to identify the velocity of the sphere, \begin{equation} U^p_i = \left( 1+\frac{a^2 }{6}\partial^2 \right)\left(u_i^\infty\left(\mathbf{x}\right)\right)| _{\mathbf{\xi}} + U^s_i. \end{equation} The angular velocity is similarly identified, \begin{equation} \Omega^p_i = \frac{1}{2}\epsilon _{ijk}\partial_j u_k^\infty| _{\mathbf{\xi}} + \Omega^s_i. \end{equation}
Supposing the position of the sphere center is given by \(X_i\) and its orientation unit vector is \(p_i\), \begin{equation} \dot{X}_i = U^p_i = \left( 1+\frac{a^2 }{6}\partial^2 \right)\left(u_i^\infty\left(\mathbf{x}\right)\right)| _{\mathbf{X}} + U^s_i, \end{equation} \begin{equation} \dot{p}_i = \epsilon _{ijk}\Omega^p_jp_k = \frac{1}{2}\epsilon _{ijk}\epsilon _{jlm}\partial_l u_m^\infty| _{\mathbf{X}}p_k + \epsilon _{ijk}\Omega^s_jp_k. \end{equation}
Method of Reflections for Widely Separated Spheres
We can solve for the flow in the presence of a single sphere; can we resolve the problem with multiple spheres? Yes, more or less. The idea behind the method of reflections is to expand the velocity field in terms of the ratio \(\alpha=a/R\), where \(a\) is the sphere’s radius and \(R\) is the separation between the spheres. We assume this ratio is small such that the expansion can be truncated at a given order and provide an accurate approximation of the multiple-sphere velocity field.
We consider two spheres centred at \(\mathbf{x}_1\) and \(\mathbf{x}_2\) with the same radius \(a\) (for simplicity; in general, the sphere radii can differ). We take the spheres to be nonrotating and translating with velocities \(U^1_i\) and \(U^2_i\). The zeroth-order solution is simply the Stokes solution for the disturbance caused by the isolated spheres, \begin{equation} u_i^{(0)}\left(\mathbf{x}\right) = u_i^{1\ (0)}\left(\mathbf{x}\right) + u_i^{2\ (0)}\left(\mathbf{x}\right) + u^\infty _i, \end{equation} and the ambient \(u^\infty _i=U^\infty_i\) is purely translation. The zeroeth-order contribution from sphere \(a\) is, \begin{equation} u_i^{a\ (0)}\left(\mathbf{x}\right) = -F_j^{a\ (0)} \left(1+\frac{a^2}{6}\partial^2\right) \frac{\mathcal{G} _{ij}\left(\mathbf{x} - \mathbf{x}^a\right)}{8\pi\mu} , \end{equation} \begin{equation} F_j^{a\ (0)} = 6\pi\mu a\left(U^\infty _j -U^a _j\right). \end{equation} We can equivalently treat the case of rotated and strained spheres using our spherical singularity expansion.
The zeroeth-order field is not quite correct since the boundary conditions at the translating spheres are not satisfied. For each translating sphere, the velocity field at the sphere surface should be the translating velocity,
\begin{equation} u_i\left(\mathbf{x} = \mathbf{x}^a\right) = U_i^a, \end{equation} In our zeroth-order expansion, the boundary condition is not satisfied due to the influence of the partner particle. For instance, at the first sphere, \begin{equation} u_i^{1\ (0)}\left(\mathbf{x} = \mathbf{x}^1\right) = U_i^1 + u_i^{2\ (0)}\left(\mathbf{x} = \mathbf{x}^1\right) + u^\infty _i. \end{equation} Our second-order corrections are the reflection fields designed to cancel the contribution from the opposing sphere. For the first particle, the first-order correction is a velocity field such that \begin{equation} u_i^{1\ (1)}\left(\mathbf{x} = \mathbf{x}^1\right) = - u_i^{2\ (0)}\left(\mathbf{x} = \mathbf{x}^1\right), \end{equation} such that the combined velocity field, \begin{equation} u_i^{(1)} = u_i^{1\ (0)} + u_i^{2\ (0)} + u_i^{1\ (1)} + u_i^{2\ (1)} \end{equation} has a boundary error that is of a higher order than the previous velocity field. This scheme can be generalized as the particle fields “reflect” on each other. For the n-th order first sphere velocity field, we find a velocity field satisfying, \begin{equation} u_i^{1\ (n)}\left(\mathbf{x} = \mathbf{x}^1\right) = - u_i^{2\ (n-1)}\left(\mathbf{x} = \mathbf{x}^1\right). \end{equation}
With our spherical singularity solution and Faxén Laws, finding reflected fields is not too difficult. For the \(a\)-th sphere, with \(b\) indices representing the opposite sphere,
\begin{equation} u_i^{a\ (1)}\left(\mathbf{x}\right) = -\left[ -F _j^{a\ (1)}\left(1+\frac{a^2}{6}\partial^2\right) + S^{a\ (1)} _{jk}\left(1+\frac{a^2}{10}\partial^2\right)\partial _k -\frac{1}{2} \epsilon _{jkl} T^{a\ (1)} _k\partial _l\right]\frac{\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}^a\right)}{8\pi\mu}. \end{equation}
The force, stress, and torque are obtained with the Faxén relations for the first sphere, \begin{equation} F_i^{1\ (1)} = 6\pi\mu a\left( 1+\frac{a^2 }{6}\partial^2 \right)u_i^{2\ (0)}| _{ \mathbf{x}=\mathbf{x}^1 }, \end{equation} \begin{equation} =-F_j^{2\ (0)}\left[ \frac{3}{4}\left(\delta _{ij}+d _id _j\right)\alpha + \frac{1}{2}\left(\delta _{ij}-3d _id _j\right)\alpha^3 \right], \end{equation} \begin{equation} T_i^{1\ (1)} = 4\pi\mu a^3 \epsilon _{ijk}\partial_ju_k^{2\ (0)}| _{ \mathbf{x}=\mathbf{x}^1 }, \end{equation} \begin{equation} = - a \alpha^2 \epsilon _{ijk}F_j^{2\ (0)}d_k, \end{equation} \begin{equation} S _{ij}^{1\ (1)} = \frac{20}{3}\pi\mu a^3\left( 1+\frac{a^2 }{10}\partial^2 \right)\frac{1}{2}\left(\partial_iu_j^{2\ (0)}+\partial_ju_i^{2\ (0)}\right)| _{ \mathbf{x}=\mathbf{x}^1 }, \end{equation} \begin{equation} = -a F_k^{2\ (0)}\left[ \frac{5}{6}\left(\delta _{ij}d_k - 3d _id _jd _k\right)\alpha^2 + \frac{4}{3}\left(-\delta _{ij}d_k-2\delta _{ik}d_j+5d _id _jd _k\right)\right]. \end{equation} Recall that \(\alpha=a/R\) is our small expansion parameter. \(R\) is the distance between the two spheres, \(R=|\mathbf{x}^1-\mathbf{x}^2|\) and \(d_i=\left(x^1_i-x_i^2\right)/R\) is a unit vector pointing between them. One can repeat this process for the second sphere and equivalently for all higher moments to expand the velocity field as a series in terms of \(\alpha\). For computations, it is often more practical to proceed in wavenumber space where one may express the reflection operation as a matrix transform.
The Effective Viscosity for a Suspension of Widely Separated Spheres
We are finally able to approach the problem of interest. The idea is as follows. Suppose we have a suspension of spherical particles in a fluid. The suspension is dilute such that the particle number density is low, and from far enough away we can treat the suspension as an ordinary fluid. This suspension fluid will not have the same viscosity as the un-suspension-ified fluid due to the particles within. Our objective is to use what we know about the interactions of particles and fluids to compute the effective viscosity of this fluid as a power series in the particle number density, which is a small number.
Formally, consider a fluid within a volume \(V\) containing a suspension of spherical particles with volume fraction \(\phi\). Let there be \(N\) particles with radius \(a\), such that the volume fraction is, \begin{equation} \phi\frac{4}{3}\pi N\frac{a^3}{V}\ll 1. \end{equation} In addition to considering a dilute suspension, we consider the limit of a widely separated suspension such that the average particle separation scale is small, \begin{equation} \langle\alpha\rangle = \frac{a}{\langle R\rangle}\ll 1. \end{equation} This follows from our assumption of a dilute suspension. As \begin{equation} \langle R \rangle \sim \left(\frac{V}{N}\right)^{\frac{1}{3}}\sim\phi^{-\frac{1}{3}}\rightarrow\alpha\sim\phi^{\frac{1}{3}}, \end{equation} and we can use our method of reflections to expand the inter-particle interactions as a power series in \(\alpha\sim\phi^{\frac{1}{3}}\).
Our objective is to compute the effective fluid viscosity \(\mu^\text{eff}\) as a power series in \(\phi\). The effective viscosity is such that the bulk average stress looks Newtonian, \begin{equation} \langle\sigma_{ij}\rangle = -\langle p\rangle \delta_{ij} + 2\mu^\text{eff}\langle e_{ij}\rangle. \end{equation} The average stress is computed by integrating over the domain \(V\), \begin{equation} \langle\sigma_{ij}\rangle = \frac{1}{V}\int_{V} dV\ \sigma_{ij} = -\langle p\rangle \delta_{ij} + 2\mu \langle e_{ij}\rangle + \frac{1}{V}\int_{\sum V_n}dV\ \sigma_{ij}, \end{equation} where we have separated the stress tensor over the fluid portion with the fluid viscosity \(\mu\) plus a contribution from averaging over all particles, where \(\sum V_n\) represents the region occupied by our suspended particles. Define this as the particle stress contribution so, \begin{equation} \sigma_{ij}^\text{particle} = \frac{1}{V}\sum_n\int_{V_n}dV\ \sigma_{ij} = \frac{1}{V}\sum_n\int_{V_n}dV\ \partial_k\left(\sigma_{ik}x_j\right) - \frac{1}{V}\sum_n\int_{V_n}dV\ x_j\partial_k\sigma_{ik} \end{equation} \begin{equation} = \frac{1}{V}\sum_n\oint_{S_n}dS\ \sigma_{ik}x_j\hat{n}_k = \frac{1}{V}\sum _n S _{ij}^{(n)} = \frac{N}{V}\langle S _{ij}\rangle. \end{equation} where \(\langle S _{ij}\rangle\) is the average stresslet over all \(N\) particles and we have used the fact that \(\partial _k\sigma _{ik}=0\) in the absence of external forcings. So all we really need to do is calculate the average stresslet, \(\langle S _{ij}\rangle\), over all particles in the suspension. We write this average as a series, \begin{equation} \langle S _{ij}\rangle = \langle S _{ij}^0\rangle + \langle S _{ij}^1\rangle + … \end{equation} Each order represents a higher order of particle-particle interactions. \(\langle S _{ij}^0 \rangle\) are the stresslet effects from isolated particles, \(\langle S _{ij}^1 \rangle\) are the stresslet effects from the first particle-particle reflection interactions, \(\langle S _{ij}^2\) are the effects from the second reflection, and so on.
Linear Corrections
At order \(\phi^0\), the particle interactions are not considered, and the stresslet is obtained from the Faxén relation, Eq. \ref{eq: stresslet flaxen relation}, \begin{equation} \langle S _{ij}^0\rangle = \frac{20}{3}\pi\mu a^3 \left(1 + \frac{a^2}{10}\partial^2\right)\langle e _{ij}\rangle. \end{equation} And to linear order, the particle stress is, \begin{equation} \frac{N}{V}\langle S _{ij}^0\rangle = \frac{20}{3}\pi\mu a^3 \frac{N}{V}\left(1 + \frac{a^2}{10}\partial^2\right)\langle e _{ij}\rangle = 5\mu\phi\langle e _{ij}\rangle + \mathcal{O}\left(\alpha^2\phi\right), \end{equation} we drop the higher-order correction since we know \(\alpha\sim\phi^{\frac{1}{3}}\), we will return to these corrections later. To linear order in \(\phi\), we obtain the effective viscosity, \begin{equation} \langle\sigma _{ij}\rangle = -\langle p\rangle \delta _{ij} + 2\mu\left(1+\frac{5}{2}\phi+\mathcal{O}\left(\phi^{5/3}\right)\right)\langle e _{ij}\rangle. \end{equation}
The Singularity System for Ellipsoids
Now we are interested in the same procedure we did for spheres but for a triaxial ellipsoid. Let the surface of the ellipsoid be given by, \begin{equation} \frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}=1, \end{equation} with \(a\geq b\geq c\). We will present a solution based on the multipole expansion, as we did for spheres.
First we define a few important quantities, \begin{equation} a_E=\sqrt{a^2-c^2}, \end{equation} \begin{equation} b_E=\sqrt{b^2-c^2}, \end{equation} \begin{equation}\label{eq: q(x) definition} q\left(\mathbf{x}\right)=\sqrt{1-\frac{x^2}{a_E^2}-\frac{y^2}{b_E^2}}. \end{equation}
Singularity System for Ellipsoids
Define the function, \begin{equation} H_n\left(\mathbf{x}\right)=\int_E dA\left(\mathbf{x’}\right)\ \frac{ q\left(\mathbf{x}’\right)^{2n-1}}{2\pi a_E b_E}\frac{1}{|\mathbf{x}-\mathbf{x}’|}. \end{equation} Because \(\partial_x^2\frac{1}{|\mathbf{x}-\mathbf{x}'|}=0\) is harmonic, \(\partial_x^2 H_n\left(\mathbf{x}\right)=0\) is also harmonic outside the region of integration by linearity. \(\int_E\) is an integral over the surface of the ellipse. Now denote the similar function, \(H_n\left(\mathbf{x};u\right)\) denoites that the particle size has been rescaled by \(u\), \(a\), \(b\), and \(c\), and replaced by \(au\), \(bu\), and \(cu\), clearly \(H_n\left(\mathbf{x}\right)=H_n\left(\mathbf{x};1\right)\). We show that \(H_{n+1}\left(\mathbf{x};1\right)\) is related to an integral over \(u\) of \(H_n\left(\mathbf{x};u\right)\). Consider, \begin{equation}\label{eq: H_n u transform} u^{2n+2}H_n\left(\mathbf{x};u\right)=\int_{E_u} dA\left(\mathbf{x’}\right)\ \frac{u^{2n}}{2\pi a_E b_E}\left(1-\frac{x^2}{u^2a_E^2}+\frac{y^2}{u^2b_E^2}\right)^{n-\frac{1}{2}}\frac{1}{|\mathbf{x}-\mathbf{x}’|}, \end{equation} \begin{equation}\label =\int_{E_u} dA\left(\mathbf{x’}\right)\ \frac{u}{2\pi a_E b_E}\left(u^2-\frac{x^2}{a_E^2}-\frac{y^2}{b_E^2}\right)^{n-\frac{1}{2}}\frac{1}{|\mathbf{x}-\mathbf{x}’|}, \end{equation} where we integrate over the scaled ellipsoid \(E_u\). The integral over \(E_u\) is over the region where, \begin{equation} \sqrt{1-\frac{x^2}{u^2a_E^2}-\frac{y^2}{u^2b_E^2}}\in\mathbb{R}, \end{equation} since \(E_u\) should be the region where the integrand is real. We can rewrite this condition as, \begin{equation} w\left(x,y\right)^2=\frac{x^2}{a_E^2}+\frac{y^2}{b_E^2}\leq u^2\leftrightarrow\left(x,y\right)\in E_u, \end{equation} which defines the region of integration \(E_u\). Now consider the integral, \begin{equation}\label{eq: w change order integration} \int_0^1 du\ \int_{E_u} dA\left(\mathbf{x’}\right)\ = \int_{E} dA\left(\mathbf{x’}\right)\ \int_w^1 du\, \end{equation} where we swapped the order of integration to integrate over \(u\) first. Performing this trick with Eq. \ref{eq: H_n u transform}, \begin{equation} \int_0^1 du\ u^{2n+2}H_n\left(\mathbf{x};u\right)=\int_0^1 du\ \int_{E_u} dA\left(\mathbf{x’}\right)\ \frac{u}{2\pi a_E b_E}\left(u^2-w^2\right)^{n-\frac{1}{2}}\frac{1}{|\mathbf{x}-\mathbf{x}’|}, \end{equation} \begin{equation} =\int _{E} dA\left(\mathbf{x’}\right)\ \int _w^1 du\ \frac{u}{2\pi a _E b _E}\left(u^2-w^2\right)^{n-\frac{1}{2}}\frac{1}{|\mathbf{x}-\mathbf{x}’|}, \end{equation}
\begin{equation} = \int_{E} dA\left(\mathbf{x’}\right)\ \frac{1}{2\pi a_E b_E}\frac{1}{2n+1}\left(1-w^2\right)^{n+\frac{1}{2}}\frac{1}{|\mathbf{x}-\mathbf{x}’|}= \int_{E} dA\left(\mathbf{x’}\right)\ \frac{1}{2\pi a_E b_E}\frac{q\left(\mathbf{x}’\right)^{2n+1}}{2n+1}\frac{1}{|\mathbf{x}-\mathbf{x}’|}, \end{equation}
so we may relate, \begin{equation}\label{eq: H_N reccurance} \left(2n+1\right)\int_0^1 du\ u^{2n+2}H_n\left(\mathbf{x};u\right)=H_{n+1}\left(\mathbf{x};1\right). \end{equation} This is neat and will be helpful.
Now consider the seemingly unrelated function, \begin{equation} G_n\left(\mathbf{x}\right)=\int_{\lambda\left(\mathbf{x}\right)}^\infty dt\ \frac{1}{\Delta\left(t\right)}p\left(x,y,z,t\right)^n, \end{equation}
where, \begin{equation} \Delta\left(t\right)=\sqrt{\left(a^2+t\right)\left(b^2+t\right)\left(c^2+t\right)}, \end{equation} \begin{equation} p\left(x,y,z,t\right)=\frac{x^2}{a^2+t}+\frac{y^2}{b^2+t}+\frac{z^2}{c^2+t}-1, \end{equation} and \(\lambda\left(\mathbf{x}\right)\) is the positive root of \(p\left(t\right)=0\). We may establish a similar relation between \(G_n\left(\mathbf{x};u\right)\) and \(G_{n+1}\left(\mathbf{x};1\right)\) as we did for \(H_n\). Consider, \begin{equation} \int_0^1 du\ u^{2n+2}G_n\left(\mathbf{x};u\right)=\int_0^1 du\ \int_{\lambda\left(\mathbf{x};u\right)}^\infty dt\ \Delta\left(t;u\right)^{-1}u^{2n+2}\left(\frac{x^2}{u^2a^2+t}+\frac{y^2}{u^2b^2+t}+\frac{z^2}{u^2c^2+t}-1\right)^n, \end{equation} now change variables with \(\tau=u^{-2} t\), and \(v\left(u\right)=\lambda\left(u\right)/u^2\), \begin{equation}\label{eq: G unchanged order integral} =\int_0^1 du\ \int_{v\left(u\right)}^\infty d\tau\ \Delta\left(\tau\right)^{-1}u\left( \frac{x^2}{a^2+\tau }+\frac{y^2}{b^2+\tau }+\frac{z^2}{ c^2+\tau }-u^2 \right)^n. \end{equation} Now we change the order of integration. We know, \begin{equation} u^2v\left(u\right)=\lambda\left(u\right)\rightarrow \frac{x^2}{u^2a^2+u^2v}+\frac{y^2}{u^2b^2+u^2v}+\frac{z^2}{u^2c^2+u^2v}=1, \end{equation} \begin{equation} \rightarrow \frac{x^2}{a^2+u^2v}+\frac{y^2}{b^2+v}+\frac{z^2}{c^2+v}=u^2. \end{equation} Therefore the inverse function of \(v\left(u\right)\), \(v^{-1}\) which is defined by, \(v^{-1}\left(v\left(u\right)\right)=u\), is given by, \begin{equation} v^{-1}\left(v\right)=\sqrt{ \frac{x^2}{a^2+v}+\frac{y^2}{b^2+v}+\frac{z^2}{c^2+v} }. \end{equation} With this, the region of integration in Eq. \ref{eq: G unchanged order integral} is, \begin{equation} 0\leq u\leq 1\quad v\left(u\right)\leq\tau\leq\infty \leftrightarrow v\left(1\right)\leq\tau\leq\infty\quad v^{-1}\left(\tau\right)\leq u \leq 1, \end{equation} and we may rewrite, \begin{equation} \int_0^1 du\ u^{2n+2}G_n\left(\mathbf{x};u\right)=\int_\lambda^\infty d\tau\ \int_{v^{-1}\left(\tau\right)}^1 du\ \Delta\left(\tau\right)^{-1}u\left( \left(v^{-1}\left(\tau\right)\right)^2-u^2 \right)^n, \end{equation} \begin{equation} =-\frac{1}{2\left(n+1\right)}\int_\lambda^\infty d\tau\ \int_{v^{-1}\left(\tau\right)}^1 du\ \Delta\left(\tau\right)^{-1} \left( \left(v^{-1}\left(\tau\right)\right)^2-u^2 \right)^{n+1}. \end{equation} so, \begin{equation}\label{eq: G_N reccurance} -2\left(n+1\right)\int_0^1 du\ u^{2n+2}G_n\left(\mathbf{x};u\right)=G_{n+1}\left(\mathbf{x};1\right). \end{equation}
Finally, we relate \(H_n\) with \(G_n\) via induction. Consider \(G_0\), \begin{equation} G_0\left(\mathbf{x}\right)=\int_{\lambda\left(\mathbf{x}\right)}^\infty dt\ \Delta\left(t\right)^{-1}. \end{equation} Now we show \(G_0\) is harmonic. \begin{equation} \partial_\lambda G_0 = -\frac{1}{\Delta\left(\lambda\right)}\rightarrow\partial_i G_0 = -\frac{1}{\Delta\left(\lambda\right)}\partial_i\lambda. \end{equation} \begin{equation} \partial^2 G_0 = -\frac{\partial^2\lambda}{\Delta\left(\lambda\right)}+\frac{\partial_\lambda\Delta\left(\lambda\right)}{\Delta\left(\lambda\right)^2}|\partial_i\lambda|^2. \end{equation} Evaluating using the definition of \(\lambda\left(\mathbf{x}\right)\), \begin{equation} \partial_i\left(\frac{x^2}{a_1^2+\lambda}+\frac{y^2}{a_2^2+\lambda}+\frac{z^2}{a_3^2+\lambda}\right)=0\rightarrow\partial_i\lambda=\frac{2x_i}{a_i^2+\lambda}\frac{1}{S\left(\lambda\right)}, \end{equation} with, \begin{equation} S\left(\lambda\right)=\frac{x^2}{\left(a^2+\lambda\right)^2}+\frac{y^2}{\left( b^2+\lambda \right)^2}+\frac{z^2}{\left(c^2+\lambda \right)^2}, \end{equation} where \(a_i=a,b,c\). This also gives us, \begin{equation} |\partial_i\lambda|^2=\frac{4}{S} \end{equation} Now for the second derivative, take another derivative, for instance along \(x\), \begin{equation} \partial_x^2\lambda=\frac{2}{a_1^2+\lambda}\frac{1}{S\left(\lambda\right)}-\frac{4x^2}{\left(a_1^2+\lambda\right)^3}\frac{1}{S\left(\lambda\right)^2}-\frac{2x}{a_1^2+\lambda}\frac{1}{S\left(\lambda\right)^2}\left( \frac{2x}{\left(a_1^2+\lambda\right)^2}-\frac{4x}{a_1^2+\lambda}\frac{T\left(\lambda\right)}{S\left(\lambda\right)}\right), \end{equation} with, \begin{equation} T\left(\lambda\right)=\frac{x^2}{\left(a^2+\lambda\right)^3}+\frac{y^2}{\left( b^2+\lambda \right)^3}+\frac{z^2}{\left(c^2+\lambda \right)^3}. \end{equation} So the laplaciacion of \(\lambda\) is, \begin{equation} \partial^2\lambda=\left(\frac{1}{a_1^2+\lambda}+\frac{1}{a_2^2+\lambda}+\frac{1}{a_3^2+\lambda}\right)\frac{2}{S\left(\lambda\right)}-\frac{4T\left(\lambda\right)}{S\left(\lambda\right)^2}-\frac{4T\left(\lambda\right)}{S\left(\lambda\right)^2}+\frac{8T\left(\lambda\right)}{S\left(\lambda\right)^2}, \end{equation} \begin{equation} =\left(\frac{1}{a_1^2+\lambda}+\frac{1}{a_2^2+\lambda}+\frac{1}{a_3^2+\lambda}\right)\frac{2}{S\left(\lambda\right)}. \end{equation} With, \begin{equation} \frac{\partial_\lambda\Delta\left(\lambda\right)}{\Delta\left(\lambda\right)}=\frac{1}{2}\left(\frac{1}{a_1^2+\lambda}+\frac{1}{a_2^2+\lambda}+\frac{1}{a_3^2+\lambda}\right), \end{equation} we see, \begin{equation} \partial^2\lambda=\frac{\partial_\lambda\Delta\left(\lambda\right)}{\Delta\left(\lambda\right)}|\partial_i\lambda|^2, \end{equation} from which we find that \begin{equation} \partial^2 G_0 = -\frac{\partial^2\lambda}{\Delta\left(\lambda\right)}+\frac{\partial_\lambda\Delta\left(\lambda\right)}{\Delta\left(\lambda\right)^2}|\partial_i\lambda|^2 = 0, \end{equation} and \(G_0\) is biharmonic as we hoped.
Because \(G_0\) is biharmonic, we may express it in the region outside of the ellipsoid as, \begin{equation} G_0 = \int_E dA\left(\mathbf{x’}\right)\ \frac{ \sigma\left(\mathbf{x}’\right)}{|\mathbf{x}-\mathbf{x}’|}, \end{equation}
where \(\sigma\left(\mathbf{x}'\right)\) is some kernel.
We will show that \(G_0=2H_0\), obtaining \(\sigma\) via Gauss’s law, \begin{equation} \left(\frac{\partial G_0}{\partial n}\right)^+-\left(\frac{\partial G_0}{\partial n}\right)^-=-4\pi\sigma, \end{equation} where \(\left(\frac{\partial G_0}{\partial n}\right)^\pm\) is the derivative of \(G_0\) in the normal direction approaching the ellipsoid surface from either the outside or the inside of the region. Now, \begin{equation} \frac{\partial G_0}{\partial n}=\frac{\partial G_0}{\partial\lambda}\frac{\partial\lambda}{\partial n}. \end{equation} Using our previous computations, \begin{equation} \hat{n}_i=\frac{\partial_i\lambda}{|\partial_i\lambda|}\rightarrow\frac{\partial\lambda}{\partial n}=\hat{n}_i\partial_i\lambda=|\partial_i\lambda|=\frac{2}{\sqrt{S\left(\lambda\right)}}, \end{equation} so, \begin{equation} \frac{\partial G_0}{\partial n}=-\frac{2}{\Delta\left(\lambda\right)}\frac{1}{\sqrt{S\left(\lambda\right)}}. \end{equation} \(\frac{\partial G_0}{\partial n}\) is discontinuous as \({\Delta\left(\lambda\right)}\rightarrow0\), we pick \(\lambda=-c^2+\epsilon\) and take the limit as \(\epsilon=0\). To leading order, \begin{equation} \frac{1}{\Delta\left(\lambda\right)} = \frac{1}{\sqrt{\left(a_E^2+\epsilon\right)\left(b_E^2+\epsilon\right)\epsilon}}\approx\frac{1}{a _Eb _E\epsilon}. \end{equation} Now, since by the definition of \(\lambda\), \begin{equation} \frac{x^2}{a^2+\lambda}+\frac{y^2}{b^2+\lambda}+\frac{z^2}{c^2+\lambda}=1\rightarrow\frac{z^2}{\epsilon} = 1 - \frac{x^2}{a_E^2} - \frac{y^2}{b_E^2} = q^2. \end{equation} So, \begin{equation} \frac{1}{\sqrt{S\left(\lambda\right)}} = \left(\frac{x^2}{\left( a_E^2+\epsilon \right)^2}+\frac{y^2}{\left( b_E^2+\epsilon \right)^2}+\frac{z^2}{\epsilon^2}\right)^{-1/2}\approx\frac{\epsilon}{q}=\frac{\sqrt{\epsilon}}{q}. \end{equation} This gives us, \begin{equation} \frac{\partial G_0}{\partial n}\approx-\frac{2}{a_Eb_Eq}, \end{equation} as \(\epsilon\rightarrow0\). Taking then the derivative from the inside of the sphere in the same way, \(\left(\frac{\partial G_0}{\partial n}\right)^+=-\left(\frac{\partial G_0}{\partial n}\right)^-\), we obtain, \begin{equation} \sigma = \frac{1}{\pi}\frac{1}{a_Eb_Eq}. \end{equation} This value of \(\sigma\) yeilds, \begin{equation} G_0 = 2\int_E dA\left(\mathbf{x’}\right)\ \frac{1}{2\pi a_E b_E}\frac{ q\left(\mathbf{x}’\right)^{-1}}{|\mathbf{x}-\mathbf{x}’|} = 2H _0. \end{equation}
With this established, all we need to do is consult our recurrence relations. From Eq. \ref{eq: H_N reccurance} and Eq. \ref{eq: H_N reccurance}, \begin{equation} H_1=\frac{1}{2}\int_0^1du\ u^{2n+2}G_0=-\frac{1}{4}G_1. \end{equation} Assume the following holds, \begin{equation}\label{eq: Hn Gn relation} H_n=\frac{\left(-1\right)^{n}}{2^{2n+1}}\frac{\left(2n\right)!}{n!n!}G_n, \end{equation} which satisfies the base case of \(n=1\). We demonstrated this holds for \(n=1\) to unite \(H_n\) and \(G_n\) with induction. Multiply by \(u^{2n+2}\) and integrate, \begin{equation} \int_0^1du\ u^{2n+2}H_n=\frac{\left(n\right)^{-1}}{2^{2n+1}}\frac{\left(2n\right)!}{n!n!}\int_0^1du\ u^{2n+2} G_n, \end{equation} \begin{equation} H_{n+1}=\frac{\left(-1\right)^{n}}{2^{2n+1}}\frac{\left(2n\right)!}{n!n!}\frac{\left(2n+1\right)}{-2\left(n+1\right)} G_{n+1}, \end{equation} \begin{equation} H_{n+1}=\frac{\left(-1\right)^{n+1}}{2^{2\left(n+1\right)}}\frac{\left(2\left(n+1\right)\right)!}{\left(n+1\right)!\left(n+1\right)!} G_{n+1}, \end{equation} proving our relation with induction.
The big key takeaway here is that since \(H_n\left(\mathbf{x}\right)\) is harmonic for all \(n\), we have proven that \(G_n\left(\mathbf{x}\right)\) is also harmonic for all \(n\). We will also use this identity to relate the \(G_n\) functions to the usual \(\mathcal{G}_{ij}\) Oseen tensor.
Multipole Expansion for Ellipsoids
We are prepared to actually build the velocity field \(u_i\). First a few derivatives of \(G_n\), \begin{equation} \partial_iG_n=-\frac{1}{\Delta\left(\lambda\right)}p\left(\lambda\right)^n\partial_i\lambda+\int_\lambda^\infty dt\ \frac{1}{\Delta\left(t\right)}\partial_i\left(p\left(t\right)^n\right), \end{equation} \begin{equation} =\int_\lambda^\infty dt\ \frac{1}{\Delta\left(t\right)}\frac{2nx_i}{a_i^2+t}p\left(t\right)^{n-1}, \end{equation} where the first term vanished since \(p\left(\lambda\right)=0\) by definition. Similarly, \begin{equation} \partial_i\partial_jG_n=2n\delta_{ij}\int_\lambda^\infty dt\ \frac{1}{\Delta\left(t\right)}\frac{1}{a_i^2+t}p\left(t\right)^{n-1}+4n\left(n-1\right)x_ix_j\int_\lambda^\infty dt\ \frac{1}{\Delta\left(t\right)}\frac{p\left(t\right)^{n-2}}{\left(a_i^2+t\right)\left(a_j^2+t\right)}, \end{equation} and for an arbitrary number of derivatives,
\begin{equation} \partial_{k_1}\partial_{k_2}…\partial_{k_n}G_n\nonumber \end{equation} \begin{equation} =\sum_{m=0}^{\lfloor n/2 \rfloor}2^{n-m}\frac{n!}{m!}\left[\left(\delta_{k_1 k_2}…\delta_{k_{2m-1}k_{2m}}x_{k_{2m+1}}…x_{k_n}\right) \int_\lambda^\infty dt\ \frac{1}{\Delta\left(t\right)}\frac{p\left(t\right)^{m} }{\prod_{\alpha=1}^n \left(a_{k_\alpha}^2+t\right)} + … + \left(\text{sym}\right) \right]. \end{equation}
In the term shown, \(\alpha\neq k_2,k_4,...,k_{2m}\). The \(\left(\text{sym}\right)\) denotes that at each \(m\) we take all permutations of the indices \(k_i\), and we have only shown one explicitly. The key idea is that at the ellipsoid surface, \(\lambda=0\), the definite integrals become constants depending on the shape of the ellipse. So at the ellipsoid surface with \(\lambda=0\), we find,
\begin{equation} \partial_{k_1}\partial_{k_2}…\partial_{k_n}\partial_{k_n+1}G_n\vert _{\lambda=0} = -2n!x _{k_1}x _{k_2}…x _{k_n}\frac{1}{\Delta\left(0\right)\prod _{\alpha=1}^n \left(a _{k _\alpha}^2\right)}\frac{\partial\lambda}{\partial x _{k _{n+1}}}\vert _{\lambda=0} \nonumber \end{equation} \begin{equation}\label{eq: partial n+1 Gn} +\text{polynomial of degree } n-1. \end{equation}
At the surface,
\begin{equation} \frac{\partial\lambda}{\partial x _{k _{n+1}}}\vert _{\lambda=0} = \frac{2x _{k _{n+1}}} { a _{k _{n+1} }^2 }\sum _{l=1}^3\left(\frac{x_l^2}{a _l^4}\right)^{-1}. \end{equation}
Thus Eq. \ref{eq: partial n+1 Gn} has no term of order \(n\), only a term of order \(n+1\) and terms of \(n-1\) and lower. We now compress the derivative notation such that \(\partial _{k_1}\partial _{k_2}...\partial _{k_n}=\partial^n\). Observe the leading-order terms at the particle surface,
\begin{equation} -x_j\partial_i\partial^{n-1}G _{n-1}| _{\lambda=0} = 2^n \left(n-1\right)! x_i x_j x _{k_1}x _{k_2}…x _{k _{n-1}}a_i^{-2}\left[\Delta\left(0\right)\left( \sum _{l=1}^3\frac{ x_l^2 }{a _l^4} \right) \prod _{\alpha=1}^{n-1} \left( a _{k _\alpha}^2 \right)\right]^{-1}+…, \end{equation}
\begin{equation} \partial_i\partial_j\partial^{n-1}G _n| _{\lambda=0} = -2^{n+1} n! x_i x_j x _{k_1}x _{k_2}…x _{k _{n-1}}a_i^{-2}a_j^{-2}\left[\Delta\left(0\right)\left( \sum _{l=1}^3\frac{x_l^2 }{a _l^4} \right) \prod _{\alpha=1}^{n-1} \left(a _{k _\alpha}^2\right)\right]^{-1}+…, \end{equation} so the combination,
\begin{equation}\label{eq: superposition G functions for v} A\delta _{ij}\partial^{n-1}G _{n-1} - x_j\partial_i\partial^{n-1} G _{n-1} + \frac{a_j^2 }{2n}\partial_i\partial_j\partial^{n-1}G _n, \end{equation}
has only terms of order \(n-1\) and lower when evaluated on the particle surface; the order \(n+1\) terms cancel entirely. The coefficient \(A\) can be picked to cancel the \(n-1\) coefficient of the second two terms while simultaneously prescribing a boundary condition at order \(n-1\).
Now for the expansion. At the particle surface, suppose an \(\left(n-1\right)\)-th order field is prescribed, so the flow at the particle is given by, \begin{equation} u_i^p=Q^{(n-1)} _{i k_1 k_2 …k _{n-1}} x _{k_1}x _{k_2}x _{k _(n-1)}. \end{equation} The velocity field requires a polynomial of degree \(n-1\) to represent the boundary condition, so it must contain a term of the form, \begin{equation} u_i \ni A^{(n-1)} _{jk _1 k_2…k _n}\partial _{k_1}\partial _{k_2}…\partial _{k _{n-1}}\left[\left(\delta _{ij}-x_j\partial_i\right)G _{n-1} + \frac{a_j^2}{2n}\partial_i\partial_jG _n\right], \end{equation} which, on the boundary, contains terms of order \(n-1\), \(n-3\), \(n-5\), and so forth. The coefficient \(\mathbf{A}^{(n-1)}\) is picked to match the boundary condition \(\mathbf{Q}^{(n)}\). To cancel the order \(n-3\) terms, the velocity field must contain a term that is of order no higher than \(n-3\) on the boundary, \begin{equation} u_i \ni A^{(n-3)} _{jk _1 k_2…k _{n-3}}\partial _{k_1}\partial _{k_2}…\partial _{k _{n-3}}\left[\left(\delta _{ij}-x_j\partial_i\right)G _{n-3} + \frac{a_j^2}{2\left(n-2\right)}\partial_i\partial_jG _{n-2}\right], \end{equation} where \(A^{(n-3)} _{jk _1 k_2...k _{n-3}}\) is picked to cancel the \(n-3\) terms from the higher order first polynomial. We repeat this process to represent the velocity field as a sum of these polynomials,
\begin{equation}\label{u solution Gn rep} u_i = \sum_{m=0}^{\lfloor \left(n-1\right)/2 \rfloor} \left(2\left(n-2m\right)-1\right) C\left(n,m\right)\frac{L _{(n-2m-1)j}}{8\pi\mu} \left[\left(\delta _{ij}-x_j\partial _i\right)G _{n-2m-1} + \frac{a_j^2}{2\left(n-2m\right)}\partial_i\partial_jG _{n-2m}\right], \end{equation} we pulled out \(\left(2\left(n-2m\right)-1\right)\) for convenience, \(C\left(n,m\right)\) is a coefficient we will discuss shortly, and the \(\mathbf{L}\) operator is defined as, \begin{equation} L _{(p)j}= \left(-1\right)^n A^{(p)} _{jk_1k_2…k _{p}}\partial _{k_1}\partial _{k_2}…\partial _{k _{p}}. \end{equation} This is excellent; we have a series solution for the velocity field provided that the \(A^{(n)} _{jk_1k_2...k _{n}}\) coefficients can be solved for. We want to derive a multipole expansion, which we can do by relating our \(G_n\) terms to \(H_n\) terms using Eq. \ref{eq: Hn Gn relation} and then relating the \(H _n\) terms to the Oseen tensor \(\mathcal{G} _{ij}\).
Combining Eq. \ref{eq: Hn Gn relation} with Eq. \ref{u solution Gn rep},
\begin{equation} u_i = \sum_{m=0}^{\lfloor \left(n-1\right)/2 \rfloor} \left(2\left(n-2m\right)-1\right)C\left(n,m\right) \left( \left(-1\right)^{n-2m-1}2^{2\left(n-2m\right)-1}\frac{\left(n-2m-1\right)!\left(n-2m-1\right)!}{\left(2\left(n-2m-1\right)\right)!}\right) \nonumber \end{equation}
\begin{equation} \times \frac{L _{(n-2m-1)j}}{8\pi\mu} \left[\left(\delta _{ij}-x_j\partial _i\right)H _{n-2m-1} - \left(2\frac{\left(n-2m\right)}{2\left(n-2m\right)-1}\right)\frac{a_j^2}{2\left(n-2m\right)}\partial_i\partial_jH _{n-2m}\right]. \end{equation}
Now, to relate to the Oseen tensor. First consider, \begin{equation} \mathcal{G}_{ij}\left(\mathbf{x}-\mathbf{x}’\right)=\frac{\delta _{ij}}{|\mathbf{x}-\mathbf{x}’|}-\left(\mathbf{x}-\mathbf{x}’\right)_i\partial _j\left(\frac{1}{|\mathbf{x}-\mathbf{x}’|}\right), \end{equation} and, \begin{equation} \partial^2\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right)=-2\partial_i\partial_j\left(\frac{1}{|\mathbf{x}-\mathbf{x}’|}\right). \end{equation} Now using the definition of \(q\left(\mathbf{x}'\right)\), Eq. \ref{eq: q(x) definition}, \begin{equation} \frac{\partial q\left(\mathbf{x}’\right)^{p+2}}{\partial x’_j}=-\frac{\left(p+2\right)x_j’}{a_j^2-c^2}q\left(\mathbf{x}’\right)^{p}, \end{equation} such that we may write, \begin{equation} \int_E dA’\ x_j’\frac{q\left(\mathbf{x}’\right)^p}{|\mathbf{x}-\mathbf{x}’|} = -\frac{a_j^2-c^2}{n+2}\int_E dA’\ \frac{1}{|\mathbf{x}-\mathbf{x}’|}\frac{\partial q\left(\mathbf{x}’\right)^{p+2}}{\partial x’_j}. \end{equation} Integrating by parts, \begin{equation} \int_E dA’\ x_j’\frac{q\left(\mathbf{x}’\right)^p}{|\mathbf{x}-\mathbf{x}’|} =\frac{a_j^2-c^2}{n+2}\int_E dA’\ q\left(\mathbf{x}’\right)^{p+2} \frac{\partial }{\partial x’_j}\left[\frac{1}{|\mathbf{x}-\mathbf{x}’|}\right]\nonumber \end{equation} \begin{equation} \int_E dA’\ x_j’\frac{q\left(\mathbf{x}’\right)^p}{|\mathbf{x}-\mathbf{x}’|} = -\frac{a_j^2-c^2}{p+2}\frac{\partial}{\partial x_j}\int_E dA’\ q\left(\mathbf{x}’\right)^{p+2} \frac{1}{|\mathbf{x}-\mathbf{x}|}, \end{equation}
where we have removed the boundary terms. With the defition of our \(H_n\) functions, we have, \begin{equation} \left(\delta _{ij}-x_j\partial _i\right)H _{n-2m-1} - \frac{a_j^2}{2\left(n-2m\right)-1}\partial_i\partial_jH _{n-2m}\nonumber \end{equation} \begin{equation} = \frac{1}{2\pi a_E b_E}\int_E dA’\ q\left(\mathbf{x}’\right)^{2\left(n-2m\right)-3}\left[\left(\delta _{ij}-x_j\partial _i\right) - \frac{a_j^2}{2\left(n-2m\right)-1}q\left(\mathbf{x}’\right)^2\partial_i\partial_j \right]\frac{1}{|\mathbf{x}-\mathbf{x}’|}. \end{equation}
Relating the first term to the Oseen tensor,
\begin{equation} \int_E dA’\ \left[q\left(\mathbf{x}’\right)^{2\left(n-2m\right)-3}\left(\delta _{ij}-x_j\partial _i\right)\right]\frac{1}{|\mathbf{x}-\mathbf{x}’|} = \int_E dA’\ q\left(\mathbf{x}’\right)^{2\left(n-2m\right)-3}\left( \mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right) - x’_j\partial_i\left[\frac{1}{|\mathbf{x}-\mathbf{x}’|} \right]\right) , \end{equation}
\begin{equation} = \int_E dA’\ q\left(\mathbf{x}’\right)^{2\left(n-2m\right)-3} \mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right) + \frac{a_j^2-c^2}{2\left(n-2m\right)-1}\partial_i\partial_j \int_E dA’\ q\left(\mathbf{x}’\right)^{2\left(n-2m\right)-1} \frac{1}{|\mathbf{x}-\mathbf{x}’|}, \end{equation}
\begin{equation} =\int_E dA’\ q\left(\mathbf{x}’\right)^{2\left(n-2m\right)-3} \mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right) - \frac{a_j^2-c^2}{4\left(n-2m\right)-2}\int_E dA’\ q\left(\mathbf{x}’\right)^{2\left(n-2m\right)-1} \partial^2\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right). \end{equation} For the second term, \begin{equation} -\int_E dA’\ \frac{a_j^2}{2\left(n-2m\right)-1}q\left(\mathbf{x}’\right)^{2\left(n-2m\right)-1}\partial_i\partial_j \left[\frac{1}{|\mathbf{x}-\mathbf{x}’|}\right]\nonumber \end{equation} \begin{equation} = \int_E dA’\ \frac{a_j^2}{4\left(n-2m\right)-2}q\left(\mathbf{x}’\right)^{2\left(n-2m\right)-1}\partial^2\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right), \end{equation} Hence, \begin{equation} \left(\delta _{ij}-x_j\partial _i\right)H _{n-2m-1} - \frac{a_j^2}{2\left(n-2m\right)+1}\partial_i\partial_jH _{n-2m} \nonumber \end{equation} \begin{equation} = \frac{1}{2\pi a_E b_E}\int_E dA’\ q\left(\mathbf{x}’\right)^{2\left(n-2m\right)-3}\left[ 1 + \frac{c^2q\left(\mathbf{x}’\right)^2}{4\left(n-2m\right)-2}\partial^2 \right] \mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right). \end{equation}
And we write down the velocity field using the Oseen tensor,
\begin{equation} u_i = \sum_{m=0}^{\lfloor\left(n-1\right)/2 \rfloor} \frac{1}{2\pi a_E b_E} L _{(n-2m-1)j} \int_E dA’\ \left(2\left(n-2m\right)-1\right) q\left(\mathbf{x}’\right)^{2\left(n-2m\right)-3}\left[ 1 + \frac{c^2q\left(\mathbf{x}’\right)^2}{4\left(n-2m\right)-2}\partial^2 \right] \frac{\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right)}{8\pi\mu}, \end{equation} where we have defined, \begin{equation} C\left(n,m\right) = \frac{\left(-1\right)^{n-2m-1}}{2^{2\left(n-2m\right)-1}}\frac{\left(2\left(n-2m-1\right)\right)!}{\left(n-2m-1\right)!\left(n-2m-1\right)!}, \end{equation} for convenience. All that remains is determing the \(\mathbf{A}^{(n)}\) coefficients to the \(\mathbf{Q}^{(n)}\) boundary condition.
Far-Field Series Form For Ellipsoidal Velocity Field
Define the following integral, \begin{equation}\label{eq: Jn integral} J_{n\ ij}\left(\mathbf{x}\right) = \int_E dA’\ f_n\left(\mathbf{x}’\right) \left[ 1 + \frac{c^2q\left(\mathbf{x}’\right)^2}{4n-2}\partial^2 \right] \frac{\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right)}{8\pi\mu}, \end{equation} with, \begin{equation} f_n\left(\mathbf{x}’\right) = \frac{2n-1}{2\pi a_E b_E} q\left(\mathbf{x}’\right)^{2n-3} = \frac{2n-1}{2\pi a_E b_E}\left(1-\frac{x’^2}{a_E^2}-\frac{y’^2}{b_E^2}\right)^{n-3/2}. \end{equation} The velocity field can be expressed as,
\begin{equation} u_i = \sum_{m=0}^{\lfloor\left(n-1\right)/2 \rfloor} \frac{L _{(n-2m-1)j}}{8\pi\mu} J _{\left(n-2m\right)\ ij}. \end{equation}
Eq. \ref{eq: Jn integral} is certainly not nice to write or study; for our purposes, it would be much more practical to obtain a series representation of this integral that we can evaluate in the far-field limit as an asymptotic solution.
First, consider, \begin{equation} I _{ij}\left(\mathbf{x}\right)=\int_E dA’\ f_1\left(\mathbf{x}’\right) \mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right). \end{equation} The integration region is where the function \(f_1\left(\mathbf{x}'\right)\) is real, \begin{equation} \left(x,y\right)\in E\leftrightarrow \frac{x^2}{a_E^2}+\frac{y^2}{b_E^2}\leq1. \end{equation} Now consider the change of variables, \begin{equation} x’=a_E\rho’\cos\phi’\quad y’=b_E\rho’\sin\phi’, \end{equation} so the integral becomes, \begin{equation} I _{ij}\left(\mathbf{x}\right)= \int_0^{2\pi}\int_0^1 d\phi’\ d\rho\ \frac{\rho}{2\pi\sqrt{1-\rho^2}}\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right). \end{equation} Now substitute \(\rho=\sin\theta\), \begin{equation} I _{ij}\left(\mathbf{x}\right)= \frac{1}{2\pi}\int_0^{2\pi}\int_0^{\pi/2} d\phi’\ d\theta’\ \sin\theta’\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right). \end{equation} Our integral has been converted to a simple integration over the upper half unit sphere. Now turning our attention to the Green’s function. Using our changed variables, \begin{equation}\label{eq: G exponent form} \mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right) = e^{-\left(x’\partial_x+y’\partial_y\right)}\mathcal{G} _{ij}\left(\mathbf{x}\right) = e^{-\sin\theta’\left(a_E\cos\phi’\partial_x+b_E\sin\phi’\right)}\mathcal{G} _{ij}\left(\mathbf{x}\right). \end{equation} Now define, \begin{equation} \tilde{D}^2=a_E^2\partial^2_x+b_E^2\partial^2_y, \end{equation} and also the angle, \begin{equation} \alpha=\tan^{-1}\left(\frac{a_E\partial_x}{b_E\partial_y}\right), \end{equation} we may transform the exponent in Eq. \ref{eq: G exponent form} to, \begin{equation} e^{-\sin\theta’\left(a_E\cos\phi’\partial_x+b_E\sin\phi’\right)} = e^{-\tilde{D}\sin\theta’\sin\left(\phi’+\alpha\right)}. \end{equation} The integral operator becomes, \begin{equation} I _{ij}\left(\mathbf{x}\right)= \frac{1}{2\pi}\int_0^{2\pi}\int_0^{\pi/2} d\phi’\ d\theta’\ \sin\theta’ e^{-\tilde{D}\sin\theta’\sin\left(\phi’+\alpha\right)} \mathcal{G} _{ij}\left(\mathbf{x}\right). \end{equation} Because we integrate from \(\phi'=0\) to \(\phi'=2\pi\) and the function \(\sin\left(\phi'+\alpha\right)\) is periodic in this range, we may remove \(\alpha\) without changing the value of the integral, \begin{equation} I _{ij}\left(\mathbf{x}\right)= \frac{1}{2\pi}\int_0^{2\pi}\int_0^{\pi/2} d\phi’\ d\theta’\ \sin\theta’ e^{-\tilde{D}\sin\theta’\sin\left(\phi’\right)} \mathcal{G} _{ij}\left(\mathbf{x}\right). \end{equation} We notice the integrand is equivalent when integrated in the range \(\pi/2\leq\theta'\leq\pi\), so we convert the integral to act over the complete unit sphere, \begin{equation} I _{ij}\left(\mathbf{x}\right)= \frac{1}{4\pi}\int_0^{2\pi}\int_0^{\pi} d\phi’\ d\theta’\ \sin\theta’ e^{-\tilde{D}\sin\theta’\sin\phi’} \mathcal{G} _{ij}\left(\mathbf{x}\right). \end{equation} Now the fun part. The integral is over the surface of the unit sphere, and we note that the exponent, when written in Cartesian form, is \(y'\). Because we are free to orient our axis as we please due to the rotational symmetry of the region of integration, we pick a coordinate system to orient the exponent along the \(z'\) direction, such that, \begin{equation} I _{ij}\left(\mathbf{x}\right)= \frac{1}{4\pi}\int_0^{2\pi}\int_0^{\pi} d\nu’\ d\eta’\ \sin\eta’ e^{-\tilde{D}\cos\eta’} \mathcal{G} _{ij}\left(\mathbf{x}\right), \end{equation} where \(\nu'\) and \(\eta'\) are our rotated spherical coordinates. Now integrate over both variables to obtain, \begin{equation}\label{eq: I sinh D verion} I _{ij}\left(\mathbf{x}\right)= \frac{\sinh\tilde{D}}{\tilde{D}}\mathcal{G} _{ij}\left(\mathbf{x}\right). \end{equation} This is interesting; we have replaced the integral over the ellipsoid with an analytic operator.
Now for a similar identity, consider \begin{equation} I’ _{ij}\left(\mathbf{x}\right)=\int_E dA’\ f_1\left(\mathbf{x}’\right)q^2\left(\mathbf{x}’\right) \partial^2\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right). \end{equation} With the same procedure performed for our previous integral, \begin{equation} I’ _{ij}\left(\mathbf{x}\right)=\frac{1}{4\pi}\int_0^{2\pi}\int_0^{\pi} d\phi’\ d\theta’\ \sin\theta’\cos^2\theta’ e^{-\tilde{D}\sin\theta’\sin\phi’} \partial^2\mathcal{G} _{ij}\left(\mathbf{x}\right). \end{equation} We perform the same rotation trick, taking \(y'\rightarrow z'\), except now, we must also recall this requires transforming \(z'\rightarrow y'\) which modifies the term in front of the exponent, \begin{equation} I’ _{ij}\left(\mathbf{x}\right)=\frac{1}{4\pi}\int_0^{2\pi}\int_0^{\pi} d\nu’\ d\eta’\ \cos^2\nu’\sin^3\eta’e^{-\tilde{D}\cos\eta’} \partial^2\mathcal{G} _{ij}\left(\mathbf{x}\right), \end{equation} \begin{equation} I’ _{ij}\left(\mathbf{x}\right)=\frac{1}{4}\int_0^{\pi} d\eta’\ \sin^3\eta’e^{-\tilde{D}\cos\eta’} \partial^2\mathcal{G} _{ij}\left(\mathbf{x}\right). \end{equation} With a change of variables and integration by parts, \begin{equation} I’ _{ij}\left(\mathbf{x}\right)=\frac{1}{4}\int _{-1}^{1} du’\ \left(1-u’^2\right) e^{-\tilde{D}u’} \partial^2\mathcal{G} _{ij}\left(\mathbf{x}\right) = -\frac{1}{2\tilde{D}}\int _{-1}^{1} du’\ u’ e^{-\tilde{D}u’} \partial^2\mathcal{G} _{ij}\left(\mathbf{x}\right), \end{equation} from which we identify, \begin{equation}\label{eq: J sinh D verion} I’ _{ij}\left(\mathbf{x}\right)= \frac{1}{2}\frac{1}{\tilde{D}}\frac{\partial}{\partial\tilde{D}}\int _{-1}^{1} du’\ e^{-\tilde{D}u’} \partial^2\mathcal{G} _{ij}\left(\mathbf{x}\right) = \frac{1}{\tilde{D}}\frac{\partial}{\partial\tilde{D}}\left(\frac{\sinh\tilde{D}}{\tilde{D}}\right)\partial^2\mathcal{G} _{ij}\left(\mathbf{x}\right). \end{equation} With Eq. \ref{eq: I sinh D verion} and Eq. \ref{eq: J sinh D verion}, \begin{equation} J _{1\ ij}\left(\mathbf{x}\right) = \int_E dA’\ f_1\left(\mathbf{x}’\right) \left[ 1 + \frac{c^2q\left(\mathbf{x}’\right)^2}{2}\partial^2 \right] \frac{\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right)}{8\pi\mu}, \end{equation} \begin{equation} J _{1\ ij}\left(\mathbf{x}\right) =\left[\frac{\sinh\tilde{D}}{\tilde{D}} + \frac{c^2}{2}\frac{1}{\tilde{D}}\frac{\partial}{\partial\tilde{D}}\left(\frac{\sinh\tilde{D}}{\tilde{D}}\right)\partial^2\right]\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}. \end{equation}
Using the power series expansion of \(\sinh\), the operator becomes, \begin{equation}\label{eq: tilde D power series operator} \frac{\sinh\tilde{D}}{\tilde{D}} + \frac{c^2}{2}\frac{1}{\tilde{D}}\frac{\partial}{\partial\tilde{D}}\left(\frac{\sinh\tilde{D}}{\tilde{D}}\right)\partial^2 = \sum _{n=0}^\infty\left( \frac{1}{\left(2n+1\right)!}\tilde{D}^{2n}+c^2\frac{n}{\left(2n+1\right)!}\tilde{D}^{2\left(n-1\right)}\partial^2 \right). \end{equation}
If we define the operator, \begin{equation} D^2=a^2\partial_x^2+b^2\partial_y^2+c^2\partial_z^2=\tilde{D}^2+c^2\partial^2, \end{equation} we may greatly simplify Eq. \ref{eq: tilde D power series operator}. Because the Oseen tensor is biharmonic, \(\partial^2\partial^2\mathbf{\mathcal{G}}=0\), \begin{equation} \tilde{D}^{2n}=\left(D^2-c^2\partial^2\right)^{n}=D^{2n}-c^2n\tilde{D}^{2\left(n-1\right)}\partial^2. \end{equation} From this, it is not difficult to obtain the rather elegant identity, \begin{equation} \frac{\sinh\tilde{D}}{\tilde{D}} + \frac{c^2}{2}\frac{1}{\tilde{D}}\frac{\partial}{\partial\tilde{D}}\left(\frac{\sinh\tilde{D}}{\tilde{D}}\right)\partial^2 = \frac{\sinh D}{D}, \end{equation} \begin{equation} J_{1\ ij}\left(\mathbf{x}\right) = \left(\frac{\sinh D}{D}\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}. \end{equation}
The last remaining step is to relate \(J_{n\ ij}\) and \(J _{\left(n-1\right)\ ij}\), which we do with the same \(u\) trick we used to derive the relations between the \(G_n\)s and \(H_n\)s. Consider, \begin{equation} \int_0^1 du\ u^{2n}J _{n\ ij}\left(\mathbf{x};u\right) =\int_0^1 du\ \int _{E_u} dA’\ u^{2n-2}\frac{2n-1}{2\pi a_E b_E}q\left(\mathbf{x}’;u\right)^{2n-3} \left[ 1 + \frac{u^2 c^2 q\left(\mathbf{x}’;u\right)^2}{4n-2}\partial^2 \right] \frac{\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right)}{8\pi\mu}, \end{equation} \begin{equation} =\int_0^1 du\ \int _{E_u} dA’\ u\frac{2n-1}{2\pi a_E b_E}\left( u^2-\frac{ x^2 }{a_E^2}-\frac{y^2}{b_E^2}\right)^{n-3/2} \left[ 1 + \frac{c^2}{4n-2}\left( u^2-\frac{ x^2 }{a_E^2}-\frac{y^2}{b_E^2}\right)\partial^2 \right] \frac{\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right)}{8\pi\mu}. \end{equation} Changing the order of integration as we did in Eq. \ref{eq: w change order integration}, \begin{equation} = \int _{E} dA’\ \int_w^1 du\ u\frac{2n-1}{2\pi a_E b_E}\left( u^2-w^2\right)^{n-3/2} \left[ 1 + \frac{c^2}{4n-2}\left( u^2-w^2\right)\partial^2 \right] \frac{\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right)}{8\pi\mu}, \end{equation} \begin{equation} = \int _{E} dA’\ \frac{1}{2\pi a_E b_E}\left( 1-w^2\right)^{n-1/2} \left[ 1 + \frac{c^2}{4n+2}\left( 1-w^2\right)\partial^2 \right] \frac{\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right)}{8\pi\mu}. \end{equation} So that we obtain the recurrence relation, \begin{equation} J _{\left(n+1\right)\ ij}\left(\mathbf{x};1\right) = \left(2n+1\right) \int _0^1 du\ u^{2n}J _{n\ ij}\left(\mathbf{x};u\right). \end{equation} We once again employ the series expansion of \(\sinh\) and \(D\left(u\right)=uD\), \begin{equation} J _{2\ ij}\left(\mathbf{x};1\right) = 3\int_0^1 du\ u^2 \sum _{k=0}^\infty \left(\frac{1}{\left(2k+1\right)!} u^{2k}D^{2k}\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu} = 3\sum _{k=0}^\infty \left(\frac{2k+2}{\left(2k+3\right)!} D^{2k}\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}, \end{equation} \begin{equation} = 3\sum _{k=1}^\infty \left(\frac{2k}{\left(2k+1\right)!} D^{2k-2}\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu} = 3\left(\frac{1}{D}\frac{\partial}{\partial D}\right)\left(\frac{\sinh D}{D}\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}. \end{equation}
Having demonstrated the base case, we verify, \begin{equation} J _{n\ ij}\left(\mathbf{x}\right) = \left(\prod _{m=1} ^{n} \left(2m-1\right)\right) \left(\frac{1}{D}\frac{\partial }{\partial D}\right)^{n-1} \left(\frac{\sinh D}{D}\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}, \end{equation} via recursion. Assume the relation holds and integrate, \begin{equation} J _{\left(n+1\right)\ ij}\left(\mathbf{x}\right) = \left(\prod _{m=1} ^{n} \left(2m-1\right)\right) (2n+1) \int_0^1 du\ u^{2n} \left(\frac{1}{D}\frac{\partial }{\partial D}\right)^{n-1} \left(\frac{\sinh D}{D}\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}, \end{equation}
\begin{equation} = \left(\prod _{m=1} ^n \left(2m-1\right)\right) (2n+1) \int_0^1 du\ \sum _{k=n-1}^\infty u^{2n+2k-2\left(n-1\right)} \frac{\left(2k\right)\left(2k-2\right)…\left(2k-2\left(n-2\right) \right) }{\left(2k+1\right)!} D^{2k-2\left(n-1\right)}\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}, \end{equation}
\begin{equation} = \left(\prod _{m=1} ^n \left(2m-1\right)\right) (2n+1) \sum _{k=n-1}^\infty \frac{\left(2k\right)\left(2k-2\right)…\left(2k-2\left(n-2\right) \right) }{\left(2k+3\right)\left(2k+1\right)!} D^{2k-2\left(n-1\right)}\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}, \end{equation}
\begin{equation} = \left(\prod _{m=1} ^{n+1} \left(2m-1\right)\right) \sum _{k=n}^\infty \frac{2k\left(2k-2\right)\left(2k-4\right)…\left(2k-2\left(n-1\right) \right) }{\left(2k+1\right)!} D^{2k-2n}\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}, \end{equation}
\begin{equation} = \left(\prod _{m=1} ^{n+1} \left(2m-1\right)\right) \left(\frac{1}{D}\frac{\partial}{\partial D}\right)^n\left(\frac{\sinh D}{D}\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}. \end{equation}
We have finished our proof. After long last, we may write the velocity field in terms of only the derivatives of the Oseen tensor, without any integrals, \begin{equation} u_i = \sum _{m=0}^{\lfloor\left( n-1 \right)/2 \rfloor} \left(\prod _{k=1} ^{n-2m} \left(2k-1\right)\right) L _{(n-2m-1)j} \left(\frac{\partial }{\partial D}\right)^{n-2m-1} \left(\frac{\sinh D}{D}\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x} \right) }{8\pi\mu}. \end{equation}
Summary of Basic Singularity Solutions for Ellipsoids
Here we summarize the velocity field solutions in the three cases of interest for our ellipsoid, those being translation, rotation, and strain. The first of these is an order-zero boundary condition, so \(n=1\), the latter two are first order, corresponding to \(n=2\).
In ellipsoidal coordinates, \begin{equation} u_i\left(\mathbf{x}\right)-U_i^\infty=\frac{1}{16\pi\mu}A_j^{(0)}\left[ \delta_{ij}G_{0}\left(\mathbf{x}\right)-x_j\partial_iG_{0}\left(\mathbf{x}\right)+\frac{a_j^2}{2}\partial_i\partial_jG_{1}\left(\mathbf{x}\right) \right], \end{equation} \begin{equation} u_i\left(\mathbf{x}\right)-\epsilon_{ijk}\Omega^\infty_jx_k-E_{ij}^\infty x_j = \frac{3}{32\pi\mu}A_{jk}^{(1)}\partial_k\left[ \delta_{ij}G_{1}\left(\mathbf{x}\right)-x_j\partial_iG_{1}\left(\mathbf{x}\right)+\frac{a_j^2}{4}\partial_i\partial_jG_{2}\left(\mathbf{x}\right) \right]. \end{equation} In terms of the integral over the Oseen tensor, \begin{equation} u_i\left(\mathbf{x}\right)-U_i^\infty=A_j^{(0)} \int_E dA’\ f_1\left(\mathbf{x}’\right)\left[ 1 + \frac{c^2q\left(\mathbf{x}’\right)^2}{2}\partial^2 \right] \frac{\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right)}{8\pi\mu}, \end{equation}
\begin{equation} u_i\left(\mathbf{x}\right)-\epsilon _{ijk}\Omega^\infty _j x _k-E _{ij}^\infty x _j = -A _{jk}^{(1)}\partial _k \int _E dA’\ f _2\left(\mathbf{x}’\right)\left[1 + \frac{c^2q\left(\mathbf{x}’\right)^2}{6}\partial^2 \right] \frac{\mathcal{G} _{ij}\left(\mathbf{x}-\mathbf{x}’\right)}{8\pi\mu}. \end{equation}
And last but not least, in terms of our multipole series, \begin{equation}\label{eq: velocity field translation ellipsoid multipole} u _i\left(\mathbf{x}\right)-U_i^\infty=A _j^{(1)} \left( \frac{\sinh D}{D}\right) \frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}, \end{equation}
\begin{equation}\label{eq: velocity field first order torque strain ellipsoid multipole} u _i\left(\mathbf{x}\right)-\epsilon _{ijk}\Omega^\infty _j x _k-E _{ij}^\infty x _j = -3A _{jk}^{(0)}\partial _k \left(\frac{1}{D}\frac{\partial }{\partial D}\right)\left(\frac{\sinh D}{D}\right) \frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}. \end{equation}
With our representation of the velocity field, we now must identify the moments, relating the \(\mathbf{A}^{(n)}\) terms to the boundary conditions on the ellipsoid. This is done by evaluating the velocity field on the surface of the ellipsoid, then reverse-engineering the appropriate moments to satisfy the boundary conditions by inverting a linear system. This is difficult to do in the general case; instead, we focus on the limit of a prolate spheroid, where the expressions reduce to much simpler forms. A prolate spheroid has \(a>b=c\), becoming slender bodies in the limit of \(c/a\ll1\).
For a prolate spheroid, the operator becomes, \begin{equation} D^2=a^2\partial_x^2+c^2\left(\partial_y^2+\partial_z^2\right)=a_E^2\partial_x^2+c^2\partial^2, \end{equation} if we define the operator, \begin{equation} D_x^2=a_E^2\partial_x^2, \end{equation} due to the biharmonicity of the Oseen tensor, we may replace, \begin{equation} D^{2k} = D_x^{2k} + k c^2D_x^{2(k-1)} \partial^2. \end{equation} Also, \begin{equation} \frac{1}{D}\frac{\partial }{\partial D} = \frac{1}{D_x}\frac{\partial }{\partial D_x}. \end{equation} The expression in the multipole expansion can be simplified, \begin{equation} \left(\frac{1}{D}\frac{\partial }{\partial D}\right)^n\left(\frac{\sinh D}{D}\right) \mathcal{G} _{ij}\left(\mathbf{x}\right) = \left( \frac{1}{D _x}\frac{\partial }{\partial D_x} \right)^n\sum _{k=0}^\infty\frac{1}{\left(2k+1\right)!}\left(D _x^{2k}+kc^2D _x^{2k-2}\partial^2\right)\mathcal{G} _{ij}\left(\mathbf{x}\right)\nonumber \end{equation}
\begin{equation} = \left(\frac{1}{D_x}\frac{\partial }{\partial D _x} \right)^n\left(\frac{\sinh D_x}{D_x}\right)\mathcal{G} _{ij}\left(\mathbf{x}\right) + \frac{c^2}{2}\left( \frac{1}{D _x}\frac{\partial }{\partial D_x} \right)^{n+1}\left(\frac{\sinh D _x}{D _x}\right)\partial^2\mathcal{G} _{ij}\left(\mathbf{x}\right). \end{equation} Employing \(D_x\) instead of \(D\) is crucial since we may recast \(D_x\) as a translation operator to simplify evaluations of the multipole moments at the spheroid boundary.
Determining Multipole Coefficients for the Prolate Spheroid
Translation
Consider the boundary condition, \begin{equation} u_i\left(r=\text{ellipsoid surface}\right)=U^p_i. \end{equation} We evaluate the velocity field in the multipole representation, Eq. \ref{eq: velocity field translation ellipsoid multipole} to determine the relationship between \(A_j^{(0)}\) and \(U^p_i\).
First, note, \begin{equation} \left(\frac{\sinh D_x}{D_x}\right) \mathcal{G} _{ij}\left(\mathbf{x}\right) = \int _0^1 dt\ \cosh\left(tD _x\right) \mathcal{G} _{ij}\left(\mathbf{x}\right), \end{equation}
\begin{equation} = \frac{1}{2}\int_0^1 dt\ \left[e^{tD _x}+e^{-tD _x}\right]\mathcal{G} _{ij}\left(\mathbf{x}\right) = \frac{1}{2}\int _0^1 dt\ \left[ \mathcal{G} _{ij}\left(\mathbf{x}+ta_E\hat{\mathbf{e}}_x\right) + \mathcal{G} _{ij}\left(\mathbf{x}-ta_E\hat{\mathbf{e}}_x\right) \right]. \end{equation}
For the derivative term, after applying the integral product rule, we see, \begin{equation} \left(\frac{1}{D_x}\frac{\partial}{\partial D_x}\right)\left(\frac{\sinh D_x}{D_x}\right) \partial^2\mathcal{G} _{ij}\left(\mathbf{x}\right) = \int _0^1 dt\ \frac{t}{D_x}\sinh\left(tD _x\right) \partial^2\mathcal{G} _{ij}\left(\mathbf{x}\right), \end{equation} \begin{equation} = \frac{1}{2}\int _0^1 dt\ \left(1-t^2\right)\cosh\left(tD _x\right) \partial^2\mathcal{G} _{ij}\left(\mathbf{x}\right) = \frac{1}{4}\int _0^1 dt\ \left(1-t^2\right)\left[ \partial^2\mathcal{G} _{ij}\left(\mathbf{x}+ta_E\hat{\mathbf{e}}_x\right) + \partial^2\mathcal{G} _{ij}\left(\mathbf{x}-ta_E\hat{\mathbf{e}}_x\right) \right]. \end{equation} We must evaluate the integrals over \(t\) at the boundary of the elipsoid. The easiest place to evaluate on the spheroid surface is at \(x=a\) and \(y=z=0\). Because the integration is over \(t\) and not over the cartesian variables, we can subtitute these values first to simplify the expressions.
Beginning with the operators acting on the Oseen tensor, \begin{equation} \mathcal{G}_{ij}\left(\mathbf{x}=\text{ellipsoid surface}\right) = \frac{\delta _{ij}}{x} + \frac{x _ix _j}{x^3}. \end{equation} For the first part,
\begin{equation} \int_0^1 dt\ \frac{1}{\sqrt{\left(a\pm a_Et\right)^2}} =\pm \frac{1}{a _E}\log\left(2\left(x\pm a_Et\right)\right)\vert _{t=0}^{t=1}, \end{equation}
\begin{equation} \int_0^1 dt\ \cosh\left(tD\right) \frac{1}{r} = \frac{1}{2a_E}\log\left(\frac{a+a_E}{a-a_E}\right). \end{equation} We obtain, \begin{equation} \left(\frac{\sinh D_x}{D_x}\right) \frac{1}{r}\vert_\left(r=\text{ellipsoid surface}\right) = \frac{1}{2ae}\log\left(\frac{1+e}{1-e}\right), \end{equation} where \(e=a_E/a\) is the unitless spheroid eccentricity.
For the second part of the Oseen tensor, only terms parallel to the spheroid axis contribute when evaluated on the surface and result in the same integral, \begin{equation} \int_0^1dt\ \frac{\left(a\pm a_Et\right)^2}{ \left(\left(x\pm a_Et\right)^2 \right)^{3/2} } = \int_0^1dt\ \frac{1}{\left(x\pm a_Et\right)}, \end{equation} \begin{equation} \left(\frac{\sinh D_x}{D_x}\right) \frac{x^2}{r^3}\vert_\left(r=\text{ellipsoid surface}\right) = \frac{1}{2ae}\log\left(\frac{1+e}{1-e}\right). \end{equation}
In summary, at the ellipsoid surface we obtain, \begin{equation} \left(\frac{\sinh D_x}{D_x}\right) \mathcal{G} _{ij}\left(\mathbf{x}=\text{ellipsoid surface}\right) = \frac{1}{2ae}\log\left(\frac{1+e}{1-e}\right)\left(\delta _{ij} + d _id _j\right), \end{equation} where \(d_i\) is the orientation direction unit vector of the spheroid along which the longer axis is placed.
Now for the Laplacian of the Oseen tensor, \begin{equation} \partial^2\mathcal{G} _{ij}\left(\mathbf{x}=\text{ellipsoid surface}\right) = 2\frac{\delta _{ij}}{x^3} - 6\frac{x _ix _j}{x^5}. \end{equation}
The integration procedure is very similar. For convenience and clarity, we will omit a few steps. \begin{equation} \left(\frac{1}{D _x}\frac{\partial}{\partial D _x}\right)\left(\frac{\sinh D _x}{D _x}\right)\frac{1}{r^3}\vert _\left(r=\text{ellipsoid surface}\right) = \frac{1}{4a^3e^3}\left(\log\left(\frac{1-e}{1+e}\right)+\frac{2e}{1-e^2}\right). \end{equation} For the second terms of the Oseen Laplacian, only components along the principal direction are non-vanishing when evaluated on the ellipsoid surface, and we obtain, \begin{equation} \left(\frac{1}{D _x}\frac{\partial}{\partial D _x}\right)\left(\frac{\sinh D _x}{D _x}\right)\frac{x^2}{r^5}\vert _\left(r=\text{ellipsoid surface}\right) = \frac{1}{4a^3e^3}\left(\log\left(\frac{1-e}{1+e}\right)+\frac{2e}{1-e^2}\right). \end{equation} All told, with both bare and Laplacian terms we have, \begin{equation} \left(\frac{\sinh D_x}{D_x}\right) \mathcal{G} _{ij}\left(\mathbf{x}=\text{ellipsoid surface}\right) + \frac{c^2}{2}\left(\frac{1}{D _x}\frac{\partial}{\partial D _x}\right)\left(\frac{\sinh D_x}{D_x}\right)\partial^2\mathcal{G} _{ij}\left(\mathbf{x}=\text{ellipsoid surface}\right) \nonumber \end{equation} \begin{equation} = \frac{1}{2ae}\log\left(\frac{1+e}{1-e}\right)\left(\delta _{ij} + d _id _j\right) + \frac{1-e^2}{8ae^3}\left(\log\left(\frac{1-e}{1+e}\right)+\frac{2e}{1-e^2}\right)\left(2\delta _{ij} - 6d _id _j\right), \end{equation}
\begin{equation} = \frac{1}{2ae^3}\left[\left(\frac{3e^2-1}{2}L+e\right)\delta _{ij}+\left(\frac{3-e^2}{2}L-3e\right)d_id_j\right], \end{equation} where we defined \(L=\log\left(\frac{1+e}{1-e}\right)\).
Matching our boundary condition, we see,
\begin{equation} A_j^{(0)}\frac{1}{16\pi\mu ae^3 }\left[\left(\frac{3e^2-1}{2}L+e\right)\delta _{ij}+\left(\frac{3-e^2}{2}L-3e\right)d _id _j\right] =U _i^P-U _i^\infty. \end{equation}
We invert to obtain, \begin{equation} A_j^{(0)} = 6\pi\mu a \frac{8e^3}{3}\left[ \left(\frac{3e^2-1}{2}L+e \right)^{-1}\delta_{jk} - \frac{\left(3-e^2 \right)L-6e}{\left(\left(1+e^2\right)L-2e\right)\left(\left( 3e^2-1 \right)L+2e\right)}d _jd _k\right]\left(U _k^P-U _k^\infty\right). \end{equation}
it is more convenient to write \(A_j^{(0)}\) in terms parallel and perpendicular to the orientation of the ellipsoid, \begin{equation} A_j^{(0)} = 6\pi\mu a \frac{8e^3}{3}\left[ \left(\left(e^2+1\right)L-2e \right)^{-1}d_jd_k + 2\left(\left(3e^2-1\right)L+2e \right)^{-1} \left(\delta _{jk}-d_jd_k\right) \right]\left(U_k^P-U_k^\infty\right). \end{equation} From this, we identify the force on the particle as, \begin{equation} u _i\left(\mathbf{x}\right)-U_i^\infty=-F_j\left[1 + \frac{1-e^2}{2}\left(\frac{1}{D_x}\frac{\partial}{\partial D_x}\right) \partial^2 \right]\left( \frac{\sinh D _x}{D _x}\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}, \end{equation} \begin{equation} F_j = 6\pi\mu a\left[ X^Ad_jd_k + Y^A\left(\delta _{jk}-d_jd_k\right)\right]\left(U_j^\infty-U_j^P\right), \end{equation} with, \begin{equation} X^A = \frac{8e^3}{3}\left(\left(e^2+1\right)L-2e \right)^{-1}, \end{equation} and, \begin{equation} Y^A = \frac{16e^3}{3}\left(\left(3e^2-1\right)L+2e \right)^{-1}. \end{equation}
Rotation
Having determined the first-order moments, we now focus on the rotational case. \begin{equation} u_i\left(r=\text{ellipsoid surface}\right)=\epsilon_{ijk}\Omega^p_jx_k. \end{equation} The procedure is similar as with translation with somewhat more involved integrals. For brevity, and because I mainly don’t want to we simply quote the results, \begin{equation} u _i\left(\mathbf{x}\right)-\epsilon _{ijk}\Omega^\infty_jx_k =\frac{3}{2}\epsilon _{jkl}T_k\partial_l\left[1 + \frac{1-e^2}{2}\left(\frac{1}{D_x}\frac{\partial}{\partial D_x}\right) \partial^2 \right]\left(\frac{1}{D_x}\frac{\partial}{\partial D_x}\right)\left( \frac{\sinh D _x}{D _x}\right)\frac{\mathcal{G} _{ij}\left(\mathbf{x}\right)}{8\pi\mu}, \end{equation} \begin{equation}\label{eq: torque prolate spheroid} T_j = 8\pi\mu a^3\left[ X^Cd_jd_k + Y^C\left(\delta _{jk}-d_jd_k\right)\right]\left(\Omega_k^\infty-\Omega_k^P\right)-8\pi\mu a^3Y^H\epsilon _{jkl}d _ld _nE _{kn}^\infty, \end{equation} with, \begin{equation} X^C = \frac{4e^3}{3}\left( 1-e^2 \right)\left(-\left(1-e^2\right)L+2e \right)^{-1}, \end{equation} and, \begin{equation}\label{eq: Yc equation} Y^C = \frac{4e^3}{3}\left( 2-e^2 \right)\left(\left(e^2+1\right)L-2e \right)^{-1}, \end{equation} and, \begin{equation}\label{eq: Yh equation} Y^H = \frac{4e^5}{3}\left(\left(e^2+1\right)L-2e \right)^{-1}. \end{equation}
There are similar expression for strain and also higher-order moments, but those are less important.
Slender Body Theory
Consider now a needle-like rigid body with a length \(2a\) much greater than its width. We will represent the body as a distribution of stokeslets, integrating over the body to show agreement with the results of our prolate spheroid full analytic result. Consider the velocity field produced by a constant line distribution of Stokeslets along the \(x\) axis, with the velocity field given by integrating these Stokeslets, \begin{equation} u _i\left(\mathbf{x}\right) = f _j\int _{-a}^a d\xi\ \frac{\mathcal{G} _{ij}\left(\mathbf{x}+\xi\hat{\mathbf{e}}_x\right)}{8\pi\mu}. \end{equation} We evaluate along the center of the rod with \(x=0\) since the surface of the rod is moving at the same velocity Defining \(\rho=\sqrt{y^2+z^2}\approx c\ll a\), we find, \begin{equation} \int _{-a}^a d\xi\ \frac{1}{\sqrt{\xi^2+\rho^2}} = 2\log\left(\frac{2}{\epsilon}\right) + \mathcal{O}\left(\epsilon^2\right), \end{equation} where \(\epsilon=c/a\) is a small expansion parameter. Similar integrals that are part of the Oseen tensor are, \begin{equation} \int _{-a}^a d\xi\ \frac{\xi^2 }{\left(\xi^2+\rho^2\right)^{3/2}} = 2\log\left(\frac{2}{\epsilon}\right) -2 + \mathcal{O}\left(\epsilon^2\right), \end{equation} \begin{equation} \int _{-a}^a d\xi\ \frac{\xi \rho}{\left(\xi^2+\rho^2\right)^{3/2}} = \mathcal{O}\left(\epsilon^2\right), \end{equation} \begin{equation} \int _{-a}^a d\xi\ \frac{\rho^2 }{\left(\xi^2+\rho^2\right)^{3/2}} = 2 + \mathcal{O}\left(\epsilon^2\right). \end{equation} So to order \(\mathcal{O}\left(\epsilon^2\right)\), the velocity field is given by, \begin{equation} u_i = \frac{1}{4\pi\mu E}f_j\left(2d _id _j + \left(\delta _{ij}-d _id _j\right)+ \mathcal{O}\left(E\right)\right), \end{equation} where \(d_i\) is the orientation direction and \(E=1/\log\left(2/\epsilon\right)\) is a small expansion parameter. The force density depends on the translation velocity in the parallel and perpendicular directions, \begin{equation} f^\parallel = 2\pi\mu E U^\parallel, \end{equation} \begin{equation} f^\perp = 4\pi\mu E U^\perp, \end{equation} and we may write the drag per unit length as, \begin{equation} f_i\left(s\right) = -\left(c^\parallel d _id _j + c^\perp\left(\delta _{ij}-d _id _j\right)\right) U_j\left(s\right), \end{equation} where \(U_i(s)\) is the local filament velocity parametrized by length \(s\) and, \begin{equation} c^\parallel = 2\pi\mu E, \end{equation} \begin{equation} c^\perp = 4\pi\mu E. \end{equation}
These are in good agreement with our results for the prolate spheroid. Recall for the translating spheroid we obtained a force density of, \begin{equation} f_j = 3\pi\mu \left[ X^Ad_jd_k + Y^A\left(\delta _{jk}-d_jd_k\right)\right]U_i. \end{equation} the limit of a thin needle is when \(\epsilon\rightarrow1\), in this case, \begin{equation} X^A = \frac{8e^3}{3}\left(\left(e^2+1\right)L-2e \right)^{-1}\rightarrow\frac{4E}{6-3E}+\mathcal{O}\left(\epsilon^2\right), \end{equation} \begin{equation} Y^A = \frac{16e^3}{3}\left(\left(3e^2-1\right)L+2e \right)^{-1}\rightarrow\frac{8E}{6+3E}+\mathcal{O}\left(\epsilon^2\right), \end{equation} which results in the same values of \(c^\parallel\) and \(c^\perp\) for thin rigid bodies.
Jeffery Orbits
We consider a prolate spheroid in a torque free field \(T_i=0\). From Eq. \ref{eq: torque prolate spheroid}, we find rotational motion follows as, \begin{equation} \left[ X^Cd_jd_k + Y^C\left(\delta _{jk}-d_jd_k\right)\right]\left(\Omega_k^\infty-\Omega_k^P\right) = Y^H\epsilon _{jkl}d _ld _nE _{kn}^\infty, \end{equation} \begin{equation} \epsilon _{ikm}x _m \left(\Omega _k^\infty-\Omega _k^P\right) = \frac{Y^H}{Y^C}\epsilon _{ijm}x _m\epsilon _{jkl}d _ld _nE _{kn}^\infty, \end{equation} \begin{equation}\label{eq: jeffry equation prolate spheroid} \epsilon _{ikm}x _m \Omega _k^P = \epsilon _{ikm}x _m \Omega _k^\infty + \frac{Y^H}{Y^C} \left(\delta _{ik} d_n d_jx_j - d_n x_k d_i \right) E _{kn}^\infty. \end{equation} We can use this equation to understand the rotation of a torque free axisymmetric particle in a shear flow given by \(u_i^\infty = \dot{\gamma}y\hat{e}^{(x)}_i\). For such a flow, the vorticity and strain are, \begin{equation} \Omega _i^\infty = \frac{1}{2}\epsilon _{ijk}\partial_ju_k^\infty = -\frac{1}{2} \dot{\gamma}\hat{e}^{(z)}_i, \end{equation} \begin{equation} E _{ij}^\infty = \frac{1}{2}\left(\partial_iu_j^\infty + \partial_ju_i^\infty\right) = \frac{1}{2}\dot{\gamma}\left(\delta _{iy}\hat{e}^{(x)}_j + \delta _{jy}\hat{e}^{(x)}_i\right). \end{equation} To understand the time evolution of the particle orientation, simply replace \(x_i\) with \(d_i\) in Eq. \ref{eq: jeffry equation prolate spheroid}, \begin{equation} \dot{d}_i = \epsilon _{ikm}d _m \Omega _k^P = \epsilon _{ikm}d _m \Omega _k^\infty + \frac{Y^H}{Y^C} \left(\delta _{ik} d_n - d_n d_k d_i \right) E _{kn}^\infty. \end{equation} Considering each part of this expression with our torque and strain, we find, \begin{equation} \epsilon _{ikm}d _m \Omega _k^\infty = -\frac{1}{2} \dot{\gamma} \left(d _x\hat{e}^{(y)}_i - d _y\hat{e}^{(x)}_i\right), \end{equation} \begin{equation} d_j E _{ij}^\infty = \frac{1}{2}\dot{\gamma}\left(d_x\hat{e}^{(y)}_i + d_y\hat{e}^{(x)}_i\right), \end{equation} \begin{equation} d_n d_k d_i E _{kn}^\infty = d_i \dot{\gamma}d _xd _y. \end{equation} So, the evolution given component by component is, \begin{equation} \dot{d}_x = \frac{1}{2} \dot{\gamma} d_y + \frac{1}{2}\dot{\gamma}\frac{Y^H}{Y^C}\left(d_y - 2 d_x^2d_y\right), \end{equation} \begin{equation} \dot{d}_y = - \frac{1}{2} \dot{\gamma} d_x + \frac{1}{2}\dot{\gamma}\frac{Y^H}{Y^C}\left(d_x - 2 d_y^2d_x\right), \end{equation} \begin{equation} \dot{d}_z = -\dot{\gamma}d_xd_yd_z. \end{equation} If we parametrize the orientation using spherical coordinates, \begin{equation} d_x = sin\theta\cos\phi\quad d_y = \sin\theta\sin\phi\quad d_z = \cos\theta, \end{equation} we obtain, \begin{equation} -\dot{d}_z = \dot{\theta}\sin\theta = \dot{\gamma}\sin^2\theta\cos\theta\sin\phi\cos\phi, \end{equation} \begin{equation} d_x\dot{d}_y - d_y \dot{d}_x = \sin^2\theta\dot{\phi} = - \frac{1}{2} \dot{\gamma} \sin^2\theta + \frac{1}{2}\dot{\gamma}\frac{Y^H}{Y^C}\sin^2\theta\left( \sin^2\phi-\cos^2\phi\right). \end{equation} Consulting the expressions for \(Y^C\) and \(Y^H\), Eq. \ref{eq: Yc equation} and Eq. \ref{eq: Yh equation}, we may define the Bretherton constant, \begin{equation} B = \frac{Y^H}{Y^C} = \frac{e^2}{2-e^2} = \frac{1-r^2}{1+r^2}, \end{equation} where \(r=c/a\) is the spheroid aspect ratio. The angular dynamics may be solved exactly with, \begin{equation} \dot{\phi} = - \frac{1}{2} \dot{\gamma} + \frac{1}{2}\dot{\gamma}\frac{1-r^2}{1+r^2}\left( 1 - 2 \cos^2\phi\right) = -\dot{\gamma}\frac{1}{1+r^2}\left( r^2\sin^2\phi + \cos^2\phi\right). \end{equation} We solve to obtain, \begin{equation} tan\left(\phi\right) = -\frac{1}{r}\tan\left(\frac{\dot{\gamma}t}{r+r^{-1}}\right). \end{equation} The period of rotation goes with, \begin{equation} T = \frac{2\pi}{\dot{\gamma}}\left(r+r^{-1}\right). \end{equation} As the rod becomes thinner and \(r\rightarrow\infty\), the period increases and the rod becomes more aligned with the direction of flow. For a spherical object, \(r=1\), the rotation is at a constant rate.
Sources
The notes here are based on Kim and Karrila’s Microhydrodynamics Principles and Selected Applications and Eric Lauga’s Fluid Dynamics of Cell Motility, both great books with rather boring titles.