Solution (source code)

= Solution

Work in the laboratory frame, with the fluid at rest at infinity. Set
$$
f(r)=\frac{3a}{4r}+\frac{a^3}{4r^3},\qquad
g(r)=\frac{3a}{4r^3}-\frac{3a^3}{4r^5},
\qquad \mathbf u=f\mathbf U+g(\mathbf U\cdot\mathbf x)\mathbf x.
$$
At $r=a$, $f=1$ and $g=0$, so the <no-slip boundary condition> is $\mathbf u=\mathbf U$ on the translating <sphere>. As $r\to\infty$, $\mathbf u=O(r^{-1})\to0$ and the gauge-fixed <pressure> tends to zero.

It is also possible to check the exterior <Stokes equations> directly. Differentiating the radial functions gives
$$
\nabla\cdot\mathbf u
=(\mathbf U\cdot\mathbf x)\left(\frac{f'}r+rg'+4g\right)=0,
$$
and
$$
\begin{aligned}
\nabla^2\mathbf u
&=\left(f''+\frac{2f'}r+2g\right)\mathbf U
+\left(g''+\frac{6g'}r\right)(\mathbf U\cdot\mathbf x)\mathbf x\\
&=\frac{3a}{2}\left(\frac{\mathbf U}{r^3}
-\frac{3(\mathbf U\cdot\mathbf x)\mathbf x}{r^5}\right)
=\frac{\nabla p}{\mu}.
\end{aligned}
$$
Thus the <Newtonian fluid stress tensor> $\sigma=-pI+\mu(\nabla\mathbf u+\nabla\mathbf u^T)$ satisfies $\nabla\cdot\sigma=0$ throughout $r>a$.

Let $\mathbf n$ point away from the common centre on each concentric <sphere>. Applying the <divergence theorem> to the intervening fluid shell, whose inner normal is $-\mathbf n$, yields
$$
\boxed{\mathbf F_{\rm drag}
=\int_{r=a}\sigma\mathbf n\,dS
=\int_{r=R}\sigma\mathbf n\,dS,\qquad R>a.}
$$
This is the <force flux at infinity in Stokes flow>. It is legitimate even though the <stress tensor> tends to zero pointwise, because the area of the large <sphere> grows as $R^2$.

For an explicit evaluation, put $A=3a/4$ and $\mathbf n=\mathbf x/R$. The leading terms are $\mathbf u=A[\mathbf U+(\mathbf U\cdot\mathbf n)\mathbf n]/R$ and $p=2\mu A(\mathbf U\cdot\mathbf n)/R^2$. Their <velocity gradient> is
$$
\partial_j u_i
=\frac{A}{R^2}\left[-U_i n_j+U_j n_i
+(\mathbf U\cdot\mathbf n)\delta_{ij}
-3(\mathbf U\cdot\mathbf n)n_i n_j\right]
+O\!\left(\frac{a^3|\mathbf U|}{R^4}\right).
$$
After symmetrizing and subtracting the <pressure> contribution,
$$
\sigma\mathbf n
=-\frac{9\mu a}{2R^2}(\mathbf U\cdot\mathbf n)\mathbf n
+O\!\left(\frac{\mu a^3|\mathbf U|}{R^4}\right).
$$
By <rotational symmetry>, $\int n_i n_j\,d\Omega=(4\pi/3)\delta_{ij}$. Hence
$$
\mathbf F_{\rm drag}
=-\frac{9\mu a}{2}\int(\mathbf U\cdot\mathbf n)\mathbf n\,d\Omega
=\boxed{-6\pi\mu a\mathbf U}.
$$
This is <Stokes's law>, with the <force> on the <sphere> opposite its <velocity>. For a component calculation with $\mathbf U=U\mathbf e_z$, the same integral is $F_z=-(9\mu aU/2)\int\cos^2\theta\,d\Omega$.