The free-space Green's function for the Stokes flow 
The Green's functions of Stokes flow
  represent solutions of the continuity
  equation $\nabla\cdot {\bf u}=0$ and the
  singularly forced Stokes equation
  $$-\nabla P+\mu
\nabla^2{\bf u}+{\bf g}\delta({\bf x-x_0})=0 \tag{*}$$
where ${\bf g}$ is an arbitrary constant,
  ${\bf x_0}$ is an arbitrary point, and
  $\delta$ is the three-dimensional
  delta function. Introducing the Green's function ${\bf G}$, we write the solution  of (*) 
  in the form
  $$u_i({\bf x})=\frac{1}{8\pi\mu}G_{ij}({\bf x,x_0})g_j$$

The following is how to find $G_{ij}$ in the free-space infinite unbounded flow from C. POZRIKIDIS's Boundary Integral and Singularity Methods for Linearized Viscous Flow
Replace the delta function with
$$
\delta(\hat{x})=-\frac{1}{4\pi}\nabla^2\left(\frac{1}{r}\right) \tag{1}
$$
where $r=|\hat{x}|$, $\hat{x}=x-x_0$. Recalling that the pressure is a harmonic function, and balancing the dimensions of the pressure term with those of the delta function in (1), we set
$$
P=-\frac{1}{4\pi}g\cdot\nabla\left(\frac{1}{r}\right) \tag{2}
$$
Substituting (1) and (2) into (*) we obtain
$$
\mu\nabla^2 u=-\frac{1}{4\pi}g\cdot(\nabla\nabla-I\nabla^2)\left(\frac{1}{r}\right)\tag{3}
$$
Next, we express the velocity in terms of a scalar function $H$ as
$$
u=\frac{1}{\mu}g\cdot(\nabla\nabla-I\nabla^2)H \tag{4}
$$
It will be noted that the continuity equation is satisfied for any choice of $H$. Substituting (4) into (3) and discarding the arbitrary constant vector g we obtain
$$
(\nabla\nabla-I\nabla^2)\left(\nabla^2H+\frac{1}{4\pi r}\right)=0\tag{5}
$$
Clearly, (5) is satisfied by $\nabla^2H=-1/(4\pi r)$. Using (1), we have $\nabla^4H=\delta(\hat{x})$. Thus 
$$
H=-\frac{r}{8\pi} \tag{6}
$$
Substituting (6) into (4) we find
$$
u_i(x)=\frac{1}{8\pi\mu}G_{ij}(\hat{x})g_j
$$
where 
$$
G_{ij}=\frac{\delta_{ij}}{r}+\frac{\hat{x}_i\hat{x}_j}{r^3}
$$


There are several places I don't understand:
  
  
*
  
*Why can we write the pressure function $P$ in (2)? What does "Recalling... we set" mean?
  
*Why can we express the velocity $u$ as (4)? How do we know it can be of this form?
  
*What's the difference between $\nabla\nabla$ and $I\nabla^2$?
  

 A: *

*(2) is obtained by taking the divergence of the whole equation, using that $u$ is divergence free, and (formally) inverting the Laplacian (by which I mean that you get an expression like $\nabla^2 P = g\cdot\nabla (\nabla^2)\frac1r$ from which you formally "cancel" out the $\nabla^2$ from both sides).
The sentence "Recalling..." is horribly written and factually wrong. In the free Stokes flow with no source term, the pressure is a harmonic function. But when driven by a source, the pressure solve a Poisson equation and is not harmonic any more. What the author meant by "balancing dimensions" is basically counting the derivatives, and is the same as the "cancelling the Laplacian" trick I gave above. 

*Expression (4) is called an "ansatz". Basically you guess that the solution can be written in that form, and see where that logically leads you. If you don't run into any contradictions, then you are gold (but note that this only gives you one of the possibly many solutions). As to how this ansatz is made: if $\mathscr{P}$ is a partial differential operator, and $\mathscr{Q}$ is another partial differential operator, and if $\mathscr{P}\mathscr{Q} = \mathscr{Q}\mathscr{P}$ (that they commute), then a solution $u$ to the equation
$$ \mathscr{P}u = \mathscr{Q}f $$
can be found if we can write $f = \mathscr{P}h$. For in this case, $u = \mathscr{Q} h$ will solve the equation. This is an analogue of the classic integrability condition for differential forms (in which case $\mathscr{P} = \partial/\partial x$ and $\mathscr{Q} = \partial/\partial y$, which commute since they are both partial derivatives).
Note: An ansatz doesn't always work. It is sort of like puzzle-solving or detective work by inference. When it works, you write it down, wait for people to ask how you saw the ansatz, and smile smugly. When it doesn't, you toss the paper into the trash bin and hope no one finds out. :-) For example, it could be that in the above example, the function $f$ is in fact not possible to be written as $\mathscr{P}h$, in which case the ansatz will fail. A trivial example for this is if $\mathscr{P} = \mathscr{Q} = d$, the exterior derivative. Let $f$ be a one-form such that $df \neq 0$. Then we violate the integrability condition and so we cannot write $f = dh$. On the other hand, the equation is trivially solvable by setting $u = f$. 

*$\nabla\nabla$ is the double gradient. $I\nabla^2$ is the Laplacian times the identity matrix. Given a function $F$, $\nabla\nabla F$ is a matrix valued function, whose components are
$$ \left( \nabla\nabla F\right)_{ij} = \frac\partial{\partial x_i}\frac{\partial}{\partial x_j} F $$
where $I\nabla^2 F$ is also a matrix valued function whose components are gived by 
$$ \left( I\nabla^2 F\right)_{ij} = (\nabla^2 F) \delta_{ij} $$
where $\delta$ is the Kronecker symbol. 
A: Responding to (3) only, 
Griffiths does a really good job explaining the nature of the laplacian in the first chapter of Introduction to Electrodynamics. 
He uses the del ($\bigtriangledown$) operator to create a few different double combinations of the Del (gradient operator):
$$\bigtriangledown \bullet (\bigtriangledown F) = \bigtriangledown^2$$
$$\bigtriangledown \times  (\bigtriangledown F)$$
$$\bigtriangledown(\bigtriangledown F)$$
Only the first one is the Laplacian in $R^3$.   The laplacian is n-dimensional though...but I've seen $\bigtriangledown^2$ used to refer to the laplacian in n dimensions, even though I think $\bigtriangledown$ is primarily an operator in $R^3$.  In Griffiths, he describes the second one as the Curl of the gradient or something...I don't remember, but it is useful for describing magnetic fields given an electric one, or something like that. 
