Chapter 9 Stokes Flow

For an incompressible fluid, the equations for mass and momentum balance are (see §2.8 and §2.9):

$\mathbf{\nabla }\cdot 𝒖=0\text{and}\rho {\displaystyle \frac{D𝒖}{Dt}}=𝐟+\mathbf{\nabla }\cdot 𝑻,$

(9.1)

where $𝒖$ is the velocity, $\rho$ is the constant density, $t$ is time, $𝑻$ is the stress tensor and $𝐟$ is a body force (i.e. $𝐟=\rho 𝐠$).

Internal viscous stresses are produced by deformation. The fluid is Newtonian if the relationship between the rate of deformation (${\gamma }_{ij}=\partial u_i/\partial x_j$) and the stress $𝑻$ is local, linear, instantaneous and isotropic. If the fluid is also incompressible, the Newtonian constitutive law is

$𝑻=-p𝑰+\mu \left(\mathbf{\nabla }𝒖+\mathbf{\nabla }{𝒖}^T\right).$

(9.2)

Substituting this into the momentum balance (9.1) gives the incompressible Navier-Stokes equations:

$\rho {\displaystyle \frac{D𝒖}{Dt}}=-\mathbf{\nabla }p+\mu {\mathbf{\nabla }}^2𝒖+𝐟$

(9.3)

To non-dimensionalise the problem we write

$𝒖=U\widetilde{𝒖},𝒙=L\widetilde{𝒙},t={\displaystyle \frac{L}{U}}\tilde{t},p=P\tilde{p},$

(9.4)

where $U$ is a characteristic velocity scale, $L$ is a characteristic lengthscale, and $P$ is a characteristic pressure scale. Substituting these into equation (9.3) when $𝐟=\mathrm{𝟎}$ and rearranging we arrive at

${\displaystyle \frac{\rho UL}{\mu }}\left({\displaystyle \frac{\partial \widetilde{𝒖}}{\partial \tilde{t}}}+\tilde{u}\cdot \tilde{\mathbf{\nabla }}\tilde{u}\right)=-{\displaystyle \frac{PL}{\mu U}}\tilde{\mathbf{\nabla }}\tilde{p}+{\tilde{\mathbf{\nabla }}}^2\widetilde{𝒖},$

(9.5)

where $\tilde{\mathbf{\nabla }}\equiv \partial /\partial {\tilde{x}}_i$. We choose the viscous scaling for the pressure $P=\mu U/L$ (we could have chosen the alternative inertial scaling $P=\rho U^2$). This leaves one non-dimensional group called the Reynolds number:

$\mathrm{Re}={\displaystyle \frac{\rho UL}{\mu }}={\displaystyle \frac{UL}{\nu }}\equiv {\displaystyle \frac{\mathrm{inertia}}{\mathrm{viscous}\mathrm{stresses}}},$

(9.6)

where $\nu =\mu /\rho$ is the kinematic viscosity ($\mu$ is denoted the dynamic viscosity). In this section we will be focusing on the limit where $\mathrm{Re}\ll 1$ such that inertia is negligible and we can ignore the left-hand side of (9.5). The mass and momentum balance, known as Stokes Equations, are then given by

$\mathbf{\nabla }\cdot 𝒖=0\text{and}\mathrm{𝟎}=-\mathbf{\nabla }p+\mu {\mathbf{\nabla }}^2𝒖+𝐟.$

(9.7)

Taking the dot product of the full momentum equation (9.1) with the velocity $𝒖$ and integrating over a volume $\mathrm{\Omega }$ gives

${\displaystyle {\int }_{\mathrm{\Omega }}}\rho 𝒖\cdot {\displaystyle \frac{\partial 𝒖}{\partial t}}d\mathrm{\Omega }+{\displaystyle {\int }_{\mathrm{\Omega }}}\rho 𝒖\cdot (𝒖\cdot \mathbf{\nabla }𝒖)d\mathrm{\Omega }={\displaystyle {\int }_{\mathrm{\Omega }}}𝒖\cdot 𝐟d\mathrm{\Omega }+{\displaystyle {\int }_{\mathrm{\Omega }}}𝒖\cdot (\mathbf{\nabla }\cdot 𝐓)d\mathrm{\Omega }.$

(9.8)

This can also be written in summation convention

${\displaystyle {\int }_{\mathrm{\Omega }}}\rho u_i{\displaystyle \frac{\partial u_i}{\partial t}}d\mathrm{\Omega }+{\displaystyle {\int }_{\mathrm{\Omega }}}\rho u_i{u}_j{\displaystyle \frac{\partial u_i}{\partial x_j}}d\mathrm{\Omega }={\displaystyle {\int }_{\mathrm{\Omega }}}{u}_i{f}_{i}d\mathrm{\Omega }+{\displaystyle {\int }_{\mathrm{\Omega }}}{u}_i{\displaystyle \frac{\partial T_{ij}}{\partial x_j}}d\mathrm{\Omega }.$

(9.9)

Rearranging, applying the divergence theorem and mass conservation, the left-hand side becomes

	$\mathrm{LHS}$	$={\displaystyle {\int }_{\mathrm{\Omega }}}\frac{1}{2}{\displaystyle \frac{\partial }{\partial t}}(\rho u_i{u}_i)-\frac{1}{2}{u}_i{u}_i{\displaystyle \frac{\partial \rho }{\partial t}}\mathrm{d}\mathrm{\Omega }+{\displaystyle {\int }_{\mathrm{\Omega }}}\frac{1}{2}{\displaystyle \frac{\partial }{\partial x_j}}\left(\rho u_i{u}_i{u}_j\right)-\frac{1}{2}{u}_i{u}_i{\displaystyle \frac{\partial }{\partial x_j}}\left(\rho u_j\right)\mathrm{d}\mathrm{\Omega }$
		$={\displaystyle \frac{d}{dt}}{\displaystyle {\int }_{\mathrm{\Omega }}}\frac{1}{2}\rho u_i{u}_{i}d\mathrm{\Omega }+{\displaystyle {\int }_{\partial \mathrm{\Omega }}}\frac{1}{2}\rho u_i{u}_i{u}_j{n}_{j}dA-{\displaystyle {\int }_{\mathrm{\Omega }}}\frac{1}{2}{u}_i{u}_i\underset{=0}{\underbrace{\left({\displaystyle \frac{\partial \rho }{\partial t}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(\rho u_j\right)\right)}}d\mathrm{\Omega }.$

Similarly, the right-hand side becomes

	$\mathrm{RHS}$	$={\displaystyle {\int }_{\mathrm{\Omega }}}{u}_i{f}_{i}d\mathrm{\Omega }+{\displaystyle {\int }_{\mathrm{\Omega }}}{\displaystyle \frac{\partial }{\partial x_j}}(u_i{T}_{ij})-T_{ij}{\displaystyle \frac{\partial u_i}{\partial x_j}}\mathrm{d}\mathrm{\Omega }$
		$={\displaystyle {\int }_{\mathrm{\Omega }}}{u}_i{f}_{i}d\mathrm{\Omega }+{\displaystyle {\int }_{\partial \mathrm{\Omega }}}{u}_i{T}_{ij}{n}_{j}dA-{\displaystyle {\int }_{\mathrm{\Omega }}}{T}_{ij}{\displaystyle \frac{\partial u_i}{\partial x_j}}d\mathrm{\Omega }.$

The mechanical energy equation is then written as

$\underset{(1)}{\underbrace{{\displaystyle \frac{d}{dt}}{\displaystyle {\int }_{\mathrm{\Omega }}}\frac{1}{2}\rho u^{2}d\mathrm{\Omega }}}+\underset{(2)}{\underbrace{{\displaystyle {\int }_{\partial \mathrm{\Omega }}}\frac{1}{2}\rho u^2𝒖\cdot 𝐧dA}}=\underset{(3)}{\underbrace{{\displaystyle {\int }_{\mathrm{\Omega }}}𝒖\cdot 𝐟d\mathrm{\Omega }}}+\underset{(4)}{\underbrace{{\displaystyle {\int }_{\partial \mathrm{\Omega }}}{u}_i{T}_{ij}{n}_{j}dA}}-\underset{(5)}{\underbrace{{\displaystyle {\int }_{\mathrm{\Omega }}}{T}_{ij}{\displaystyle \frac{\partial u_i}{\partial x_j}}d\mathrm{\Omega }}},$

(9.10)

where

1. 1.

   the rate of change of kinetic energy in volume $\mathrm{\Omega }$ is due to:

2. 2.

   the flux of kinetic energy over the boundary $A$,

3. 3.

   the rate of working by body forces $𝐟$ in volume $\mathrm{\Omega }$,

4. 4.

   the rate of working of surface forces in $A$,

5. 5.

   and the stress working in volume $\mathrm{\Omega }$.

The mechanical energy equation above is for a general continuum. We can now look more specifically for a Newtonian viscous fluid with constitutive law,

$T_{ij}=-p{\delta }_{ij}+2\mu e_{ij}\mathrm{where}{e}_{ij}=\frac{1}{2}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right).$

(9.11)

The stress working per unit volume then becomes

$T_{ij}{\displaystyle \frac{\partial u_i}{\partial x_j}}=-p{\delta }_{ij}{\displaystyle \frac{\partial u_i}{\partial x_j}}+2\mu e_{ij}{\displaystyle \frac{\partial u_i}{\partial x_j}}.$

(9.12)

Remembering some properties of the rate of strain and rotation tensors, we know that

${\delta }_{ij}{\displaystyle \frac{\partial u_i}{\partial x_j}}={\displaystyle \frac{\partial u_i}{\partial x_i}}=\mathbf{\nabla }\cdot 𝒖=0,e_{ij}+r_{ij}=\frac{1}{2}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right)+\frac{1}{2}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}-{\displaystyle \frac{\partial u_j}{\partial x_i}}\right)={\displaystyle \frac{\partial u_i}{\partial x_j}},$

(9.13)

$e_{ij}{r}_{ij}=\frac{1}{4}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right)\left({\displaystyle \frac{\partial u_i}{\partial x_j}}-{\displaystyle \frac{\partial u_j}{\partial x_i}}\right)=0.$

(9.14)

Hence, the stress working, also known as the viscous dissipation $\mathrm{\Phi }$, per unit volume is given by

$\mathrm{\Phi }=T_{ij}{\displaystyle \frac{\partial u_i}{\partial x_j}}=2\mu e_{ij}{e}_{ij}.$

(9.15)

The viscous dissipation $\mathrm{\Phi }$ describes the rate of loss of mechanical energy due to friction in the volume $\mathrm{\Omega }$.

The approximation of Stokes Flow is good if the velocity or length scales $U,L$ are small (e.g. swimming microorganisms), or the dynamic viscosity $\mu$ is large (e.g. in geophysics, flow in porous media, or the flow of lava or ice). Some examples are:

	$L$	$U$	$\rho$	$\mu$	$\nu$	$\mathrm{Re}$
Sperm in water:	$50\mu m$	$200\mu m/s$	${10}^{3}kg/m^3$	${10}^{-3}Pas$	${10}^6{m}^2/s$	${\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$
Mantle:	$3\times {10}^{3}km$	$1cm/yr$	$3300kg/m^3$	${10}^{21}Pas$	$3\times {10}^{17}{m}^2/s$	$\mathrm{𝟑}\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐𝟏}}$

(See E. M. Purcell (1977) Life at Low Reynolds Number, American Journal of Physics).

When approximating the Reynolds number, there are some caveats. Firstly there may be more than one relevant lengthscale, e.g in Lubrication theory later in the course the vertical lengthscale $H$ is much smaller than the horizontal lengthscale $L$. These lengthscales may also vary with position e.g. flow past a body where in the far field $\mathrm{Re}\sim Ur/\nu$ where $r$ is the radial distance away from the body. Finally, the time scale may not be $L/U$ but instead may be specified by an external driving force e.g. an oscillating boundary condition with timescale $T\sim {\omega }^{-1}$.

Consider Stokes equations with body force $𝐟=\mathrm{𝟎}$:

$-\mathbf{\nabla }p+\mu {\mathbf{\nabla }}^2𝒖=\mathrm{𝟎},\mathbf{\nabla }\cdot 𝒖=0.$

(9.39)

	$1.$	$\mathbf{\nabla }\cdot \left(-\mathbf{\nabla }p+\mu {\mathbf{\nabla }}^2𝒖\right)=0\Rightarrow {\mathbf{\nabla }}^{2}p=0$		(9.40)
	$2.$	$\mathbf{\nabla }\times \left(-\mathbf{\nabla }p+\mu {\mathbf{\nabla }}^2𝒖\right)=\mathrm{𝟎}\Rightarrow {\mathbf{\nabla }}^2𝝎=\mathrm{𝟎}\mathrm{where}𝝎=\mathbf{\nabla }\times 𝒖$		(9.41)
	$3.$	${\mathbf{\nabla }}^2\left(-\mathbf{\nabla }p+\mu {\mathbf{\nabla }}^2𝒖\right)=\mathrm{𝟎}\Rightarrow {\mathbf{\nabla }}^4𝒖=\mathrm{𝟎}.$		(9.42)

Hence, $p$ and $𝝎$ are harmonic and $𝒖$ is biharmonic.

A method for finding the solution of Stokes flows comes from the Papkovich-Neuber representation. Suppose we can write the pressure $p={\mathbf{\nabla }}^2\mathrm{\Pi }$. Then

${\mathbf{\nabla }}^2(-\mathbf{\nabla }\mathrm{\Pi }+\mu 𝒖)=\mathrm{𝟎}\Rightarrow \mu 𝒖=\mathbf{\nabla }\mathrm{\Pi }-𝚽\mathrm{where}{\mathbf{\nabla }}^2𝚽=\mathrm{𝟎}.$

(9.43)

Using the incompressibility condition, we have

$\mu \mathbf{\nabla }\cdot 𝒖=\mathbf{\nabla }\cdot (\mathbf{\nabla }\mathrm{\Pi }-𝚽)=0\Rightarrow {\mathbf{\nabla }}^2\mathrm{\Pi }=\mathbf{\nabla }\cdot 𝚽\Rightarrow \mathrm{\Pi }=\frac{1}{2}(𝒙\cdot 𝚽+\chi )\mathrm{where}{\mathbf{\nabla }}^2\chi =0.$

(9.44)

We can check this last relation:

	${\displaystyle \frac{{\partial }^2\mathrm{\Pi }}{\partial x_i\partial x_i}}$	$={\displaystyle \frac{{\partial }^2}{\partial x_i\partial x_i}}\left({\displaystyle \frac{1}{2}}{x}_j{\mathrm{\Phi }}_j+{\displaystyle \frac{1}{2}}\chi \right)={\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial }{\partial x_i}}\left({\displaystyle \frac{\partial x_j}{\partial x_i}}{\mathrm{\Phi }}_j+x_j{\displaystyle \frac{\partial {\mathrm{\Phi }}_j}{\partial x_i}}\right)+{\displaystyle \frac{1}{2}}\underset{=0}{\underbrace{{\displaystyle \frac{{\partial }^2\chi }{\partial x_i\partial x_i}}}}$		(9.45)
		$={\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial }{\partial x_i}}\left({\delta }_{ij}{\mathrm{\Phi }}_j+x_j{\displaystyle \frac{\partial {\mathrm{\Phi }}_j}{\partial x_i}}\right)={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial {\mathrm{\Phi }}_i}{\partial x_i}}+{\delta }_{ij}{\displaystyle \frac{\partial {\mathrm{\Phi }}_j}{\partial x_i}}+x_j\underset{=0}{\underbrace{{\displaystyle \frac{{\partial }^2{\mathrm{\Phi }}_j}{\partial x_i\partial x_i}}}}\right)={\displaystyle \frac{\partial {\mathrm{\Phi }}_i}{\partial x_i}}.$		(9.46)

 Any Stokes flow with $𝐟=\mathrm{𝟎}$ can be written in terms of a harmonic vector $𝚽$ and a harmonic scalar $\chi$:

$2\mu 𝒖=\mathbf{\nabla }\left(𝒙\cdot 𝚽+\chi \right)-2𝚽\mathrm{and}p=\mathbf{\nabla }\cdot 𝚽.$

(9.47)

 

Spherical harmonic functions are used for solutions involving points, spheres and cylinders that have no preferred direction or orientation. Let the radial coordinate $r=|𝒙|\equiv {(𝒙\cdot 𝒙)}^{1/2}$. We can show that ${\mathbf{\nabla }}^2\frac{1}{r}=0$, except when $r=0$:

	${\displaystyle \frac{{\partial }^2}{\partial x_j\partial x_j}}\left({(x_i{x}_i)}^{-1/2}\right)$	$=-{\displaystyle \frac{\partial }{\partial x_j}}\left({(x_k{x}_k)}^{-3/2}{x}_i{\delta }_{ij}\right)=-{\displaystyle \frac{\partial }{\partial x_j}}\left({(x_k{x}_k)}^{-3/2}{x}_j\right)$		(9.48)
		$=-3{(x_k{x}_k)}^{-3/2}+3{(x_m{x}_m)}^{-5/2}{x}_k{x}_j{\delta }_{kj}=0.$		(9.49)

All other harmonic functions $\varphi$ with $\varphi \to 0$ and $r\to \mathrm{\infty }$ are made up of

${\displaystyle \frac{1}{r}},\mathbf{\nabla }{\displaystyle \frac{1}{r}},\mathbf{\nabla }\mathbf{\nabla }{\displaystyle \frac{1}{r}},\mathbf{\nabla }\mathbf{\nabla }\mathbf{\nabla }{\displaystyle \frac{1}{r}},\mathrm{\dots }$

(9.50)

Harmonic functions that are bounded at $r=0$ are made up of

$r{\displaystyle \frac{1}{r}},r^3\mathbf{\nabla }{\displaystyle \frac{1}{r}},r^5\mathbf{\nabla }\mathbf{\nabla }{\displaystyle \frac{1}{r}},\mathrm{\dots }{r}^{2n+1}{\mathbf{\nabla }}^n{\displaystyle \frac{1}{r}},\mathrm{\dots }$

(9.51)

These solutions only depend on $r$ and $𝒙$ and have no preferred orientation. Some useful derivatives to remember are:

	$1.$	${(\mathbf{\nabla }𝒙)}_{ij}={\displaystyle \frac{\partial x_i}{\partial x_j}}={\delta }_{ij}(\mathrm{identity})$		(9.52)
	$2.$	${(\mathbf{\nabla }r)}_i={\displaystyle \frac{x_i}{r}}$		(9.53)
	$3.$	${\left(\mathbf{\nabla }{\displaystyle \frac{1}{r}}\right)}_i=-{\displaystyle \frac{x_i}{r^3}}$		(9.54)
	$4.$	${(\mathbf{\nabla }f(r))}_i=f^{\prime }(r){(\mathbf{\nabla }r)}_i=f^{\prime }(r){\displaystyle \frac{x_i}{r}}$		(9.55)

Since we know that Stokes flows are linear in boundary conditions and forcing, we can form Papkovich-Neuber solutions by multiplying them by constant scalars, vectors, tensors, and forming dot-products, for example

$𝑨{\displaystyle \frac{1}{r}},B\mathbf{\nabla }{\displaystyle \frac{1}{r}},𝑪\cdot \mathbf{\nabla }{\displaystyle \frac{1}{r}},\mathrm{\dots }$

(9.56)

When building these solutions there are some properties it is useful to recall about tensors. An $n$th-rank tensor is a set of objects that transform like a product of $n$ vectors. For example, if $T_{i_1{i}_2\mathrm{\dots }{i}_n}$ is an $n$th-rank tensor, then

$T_{j_1{j}_2\mathrm{\dots }{j}_n}^{\prime }=O_{j_1{i}_1}{O}_{j_2{i}_2}\mathrm{\dots }{O}_{j_n{i}_n}{T}_{i_1{i}_2\mathrm{\dots }{i}_n},$

(9.57)

where $O_{ij}$ is a transformation matrix. Symmetry properties are tensor invariant, so if $T_{i_1{i}_2\mathrm{\dots }{i}_n}$ is symmetric (antisymmetric) under an exchange of $i_1$ and $i_2$ then $T_{j_1{j}_2\mathrm{\dots }{j}_n}^{\prime }$ is also symmetric (antisymmetric) under an exchange of $j_1$ and $j_2$. Examples of tensors in our problems would be

$\mathrm{velocity}𝒖,\mathrm{force}𝑭={\displaystyle \int 𝑻\cdot 𝐧dA},\mathrm{position}𝒙,\mathrm{grad}\mathbf{\nabla },\mathrm{identity}𝑰.$

(9.58)

In contrast, pseudotensors, such as the Levi-Civita ${ϵ}_{ijk}$, change sign under transformation by reflections. Examples of pseudotensors we will come across are

$\mathrm{angular}\mathrm{velocity}𝛀,\mathrm{couple}𝑮={\displaystyle \int 𝒙\times (𝑻\cdot 𝐧)dA},\mathrm{cross}\mathrm{product}𝒖\times 𝒙,\mathrm{vorticity}𝝎=\mathbf{\nabla }\times 𝒖.$

(9.59)

If we want to find the velocity $𝒖$ using the Papkovich-Neuber representation, then we want $𝚽$ and $\chi$ to be tensors.

Point force (Stokeslet): We want to find the velocity field for a point force. Consider the following:

	$\mathbf{\nabla }\cdot 𝑻$	$=\mu {\mathbf{\nabla }}^2𝒖-\mathbf{\nabla }p=-𝑭\delta (𝒙)$		(9.60)
	$\mathbf{\nabla }\cdot 𝒖$	$=0,𝒖\to \mathrm{𝟎}\mathrm{as}|𝒙|\to \mathrm{\infty }.$		(9.61)

We know that

1. $-$

   by linearity, the velocity field must be linear in the forcing $𝑭$

2. $-$

   solution must decay in the far field

3. $-$

   there is no preferred direction in the problem, can use spherical harmonics

4. $-$

   the velocity field $𝒖$ must be made up of tensors

Using this information, we can now try to build harmonic functions $𝚽$ and $\chi$ using the solutions in (9.50)-(9.51). Let’s try the following ansatzes:

	$𝚽$	$={\displaystyle \frac{\alpha 𝑭}{r}}\mathrm{✓}$		(9.62)
	$𝚽$	$=\alpha 𝑭\cdot \mathbf{\nabla }\mathbf{\nabla }{\displaystyle \frac{1}{r}}\left(\mathrm{equivalently}\chi =\beta 𝑭\cdot \mathbf{\nabla }{\displaystyle \frac{1}{r}}\right)\to \mathrm{too}\mathrm{singular}\mathrm{at}r=0\times$		(9.63)
	$𝚽$	$=\alpha 𝑭\times \mathbf{\nabla }{\displaystyle \frac{1}{r}}\to \mathrm{pseudotensor}\times$		(9.64)

The first ansatz is the only one that decays in the far field and will not blow up at the origin. We proceed by substituting $𝚽$ into the solution (9.47):

$2\mu 𝒖=\mathbf{\nabla }\left(\alpha {\displaystyle \frac{𝑭\cdot 𝒙}{r}}\right)-{\displaystyle \frac{2\alpha 𝑭}{r}}=-\alpha \left({\displaystyle \frac{𝑭}{r}}+{\displaystyle \frac{(𝑭\cdot 𝒙)𝒙}{r^3}}\right).$

(9.65)

To fully specify the solution we need to find $\alpha$. Going back to the problem statement, we want the solution to satisfy

$\mathbf{\nabla }\cdot 𝒖=0\Rightarrow {\displaystyle {\int }_{r=R}}𝒖\cdot 𝒏dA=0,$

(9.66)

(this is built into the Papkovich-Neuber representation), and

$-𝑭\delta (𝒙)=\mathbf{\nabla }\cdot 𝑻\Rightarrow -𝑭={\displaystyle {\int }_{r=R}}𝑻\cdot 𝒏dA.$

(9.67)

On a sphere $r=R$, $𝒙=R𝒏$, the normal velocity can be calculated

${𝒖\cdot 𝒏|}_{r=R}=-{\displaystyle \frac{\alpha }{2\mu }}\left({\displaystyle \frac{𝑭}{R}}+{\displaystyle \frac{(𝑭\cdot R𝒏)R𝒏}{R^3}}\right)\cdot 𝒏=-{\displaystyle \frac{\alpha (𝑭\cdot 𝒏)}{\mu R}}$

(9.68)

$\Rightarrow {\displaystyle {\int }_{r=R}}𝒖\cdot 𝒏dA=-{\displaystyle \frac{\alpha 𝑭}{\mu R}}\cdot \underset{=0,\mathrm{isotropic}}{\underbrace{{\displaystyle {\int }_{r=R}}𝒏dA}}=0.$

(9.69)

The normal stress can also be calculated, but requires a bit more algebra:

${𝑻\cdot 𝒏|}_{r=R}={(-p𝑰+\mu (\mathbf{\nabla }𝒖+\mathbf{\nabla }{𝒖}^T))\cdot 𝒏|}_{r=R}$

(9.70)

with

$p=\mathbf{\nabla }\cdot \left({\displaystyle \frac{\alpha 𝑭}{r}}\right)=-{\displaystyle \frac{\alpha 𝑭\cdot 𝒙}{r^3}},\mathbf{\nabla }𝒖=-{\displaystyle \frac{\alpha }{2\mu }}\left({\displaystyle \frac{(𝑭\cdot 𝒙)𝑰}{r^3}}-3{\displaystyle \frac{(𝑭\cdot 𝒙)𝒙𝒙}{r^5}}\right),$

(9.71)

${𝑻\cdot 𝒏|}_{r=R}={\displaystyle \frac{3\alpha (𝑭\cdot 𝒏)𝒏}{R^2}}\Rightarrow {{\displaystyle {\int }_{r=R}}𝑻\cdot 𝒏|}_{r=R}\mathrm{d}A={\displaystyle \frac{3\alpha 𝑭}{R^2}}\cdot \underset{𝑨}{\underbrace{{\displaystyle {\int }_{r=R}}𝒏𝒏dA}}={\displaystyle \frac{3\alpha 𝑭}{R^2}}\cdot {\displaystyle \frac{4\pi R^2𝑰}{3}}$

(9.72)

Here, we have noticed that the integral marked $𝑨$ is an isotropic second rank tensor so must be $\propto {\delta }_{ij}$. We can determine the constant of proportionality by calculating the trace, $A_{ii}=4\pi R^2$. Hence, $\alpha =-1/4\pi$. Notice that $\alpha$ is independent of $R$. If we had chosen one of the other harmonic solutions with higher powers of $1/r$, they would be too singular to satisfy the force condition.

We know that

1. $-$

   by linearity, the velocity field must be linear in the source flux $Q$

2. $-$

   solution must decay in the far field

3. $-$

   there is no preferred direction in the problem, can use spherical harmonics

4. $-$

   the velocity field $𝒖$ must be made up of tensors

Let’s try harmonic function

$\chi ={\displaystyle \frac{\beta Q}{r}}\left(\mathrm{equivalently}𝚽=\widehat{\beta }Q\mathbf{\nabla }{\displaystyle \frac{1}{r}}\right).$

(9.75)

In a similar manner to the point force, the other harmonic solutions with higher powers of $1/r$ would be too singular to satisfy the source condition. Substituting into the solution (9.47) we find:

$2\mu 𝒖=\mathbf{\nabla }\left({\displaystyle \frac{\beta Q}{r}}\right)\Rightarrow 𝒖=-{\displaystyle \frac{\beta Q𝒙}{2\mu r^3}}.$

(9.76)

To find $\beta$, we use the source condition at $r=0$:

$\mathbf{\nabla }\cdot 𝒖$

$=Q\delta (𝒙)\Rightarrow Q={\displaystyle {\int }_{r=R}}𝒖\cdot 𝒏dA=-{\displaystyle \frac{\beta Q}{2\mu R^2}}{\displaystyle {\int }_{r=R}}𝒏\cdot 𝒏dA=-{\displaystyle \frac{\beta Q}{2\mu R^2}}4\pi R^2.$

(9.77)

Hence, $\beta =-\mu /2\pi$. The velocity field then becomes $𝒖=Q𝒙/4\pi r^3$ which is independent of viscosity $\mu$. This is the three-dimensional extension of the potential flow for a point source studied in §6.

A sphere of radius $r=a$ translating with velocity $𝐔$ and rotating with angular velocity $𝛀$ moves with velocity

${𝒖|}_{r=a}=𝐔+𝛀\times 𝒙$

(9.78)

and exerts a force $𝐅$ and couple $𝐆$ on the fluid. The dissipation in the fluid due to the particle is (using work balances dissipation)

	${\displaystyle {\int }_{\mathrm{\Omega }}}\mathrm{\Phi }d\mathrm{\Omega }$	$={\displaystyle {\int }_{\partial \mathrm{\Omega }}}𝒖\cdot 𝑻\cdot 𝐧dA={\displaystyle {\int }_{particle}}(𝐔+𝛀\times 𝒙)\cdot 𝑻\cdot 𝐧dA$		(9.79)
		$=𝐔\cdot {\displaystyle {\int }_{particle}}𝑻\cdot 𝐧dA+𝛀\cdot {\displaystyle {\int }_{particle}}𝒙\times (𝑻\cdot 𝐧)dA$		(9.80)
		$=𝐔\cdot 𝐅+𝛀\cdot 𝐆,$		(9.81)

where the normal points out of the fluid into the body, with

$\mathrm{force}𝑭={\displaystyle {\int }_{particle}}𝑻\cdot 𝐧dA\mathrm{and}\mathrm{couple}𝑮={\displaystyle {\int }_{particle}}𝒙\times (𝑻\cdot 𝐧)dA$

(9.82)

exerted on the fluid by the sphere.

Velocity Field: We would like to calculate the velocity field for a rigid sphere of radius $r=a$ translating with velocity $𝒖=𝑼$. Consider the following:

	$\mathbf{\nabla }\cdot 𝑻=\mu {\mathbf{\nabla }}^2𝒖-\mathbf{\nabla }p=\mathrm{𝟎},\mathbf{\nabla }\cdot 𝒖=0,$		(9.83)
	$𝒖=𝑼\mathrm{at}r=a,𝒖\to \mathrm{𝟎}\mathrm{as}|𝒙|\to \mathrm{\infty }.$		(9.84)