The following is how to find $G_{i j}$ 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) $$
I don't understand why and how is it possible to replace the delta function in this way. Or even if there's a physical intuition (which isn't clear either) how can this match up in the distributional sense.
One finds easily that $\nabla\frac{1}{r} = -\frac{\hat{x}}{r^3}$. Calculating the divergence of this results in $\nabla\cdot\nabla\frac{1}{r} = 0$ but it should be remembered that this is only defined for $r \neq 0$. Thus, the surface integral $\oint_{\partial\Omega} \nabla\frac{1}{r} \cdot n \, dS = 0$ for any region $\Omega$ not containing origin. If one calculates $\oint_{\partial\Omega} \nabla\frac{1}{r} \cdot n \, dS = 0$ over the surface of a region containing origin, for example of a sphere with center in origin, then one gets $-4\pi$.
Thus, $\nabla \cdot \nabla\frac{1}{r} = -4\pi \, \delta$.