Examples of Greens functions for Laplace's equation with Neumann boundary conditions. The von Neumann boundary problem is a PDE in $\Omega$
\begin{cases}\Delta u＝0\\\frac{\partial u}{\partial \vec n}\rvert_{\Gamma}＝g\end{cases}where $g$ is assumed to be smooth on $\Gamma＝\partial\Omega$ $\vec n$ is the unit normal vector. $\Omega$ is bounded. 
The green function of this problem is $G(p,q)$ such that\begin{cases}   \Delta_p G(p,q)＝\delta(p-q)  \\   \frac{\partial_p G}{\partial \vec n}＝0\end{cases}
I want an explicit formula as an example. when the domain$\Omega$ is special, such as a sphere, a cube, etc.
----------05/10----------
According to Kenny Wong's answer, the conditions of $g$ should satisfy $\int_\Gamma gdS＝0$, and the boundary condition for $G(p,q)$ should be changed to $\frac{\partial G}{\partial\vec n}\lvert_\Gamma＝constant$ such.that the integral of the constant is one on $\Gamma$.
 A: I'm afraid that what I'm about to say is not quite you what to hear: your Neumann PDE does not have a solution for arbitrary choices of $g$.
To see why, let's compute the integral of $g$ over the surface $\Gamma$:
$$ \int_\Gamma g \ dS = \int_\Gamma \frac{\partial u}{\partial n} \ dS=\int_\Omega \nabla^2 u \ dV =\int_\Omega 0\ dV = 0.$$
Therefore, the average value of $g$ on the boundary surface $S$ must be zero. If $g$ fails to satisfy this condition, your PDE cannot be solved.
For a similar reason, your definition of the Green's function must  be modified. Indeed, if we integrate $\frac{\partial G}{\partial n}$ over the surface $\Gamma$, we find that
$$ \int_\Gamma \frac{\partial G(\vec r, \vec r')}{\partial n_{\vec r}}dS(\vec r) = \int_\Omega \nabla_{\vec r}^2 G(\vec r,\vec r') \ dV(\vec r) = \int_V \delta(\vec r-\vec r') \ dV(\vec r) = 1.$$
So it is inconsistent to set $\frac{\partial G}{\partial n}$ equal to zero on the boundary $\Gamma$.
Fortunately, all is not lost! We can redefine the Green's function $G$ so that it satisfies
$$ \begin{cases} \nabla_{\vec r}^2 G(\vec r,\vec r') = \delta(\vec r - \vec r') &{\rm \ \ on \ \ } \Omega , 
\\
\frac{\partial G(\vec r, \vec r')}{\partial n_{\vec r}}=\frac 1 A & {\rm \ \ on \ \ } \Gamma, \end{cases}
$$
where $A = \int_\Gamma dS$ is the area of the boundary surface $\Gamma$.
Now, Green's identity states that
\begin{multline}\int_\Omega \left( u(\vec r) \nabla_{\vec r}^2 G(\vec r , \vec r') - G(\vec r,\vec r') \nabla_{\vec r} u(\vec r) \right) \ dV(\vec r) \\ = \int_\Gamma \left( u(\vec r) \frac{\partial G(\vec r,\vec r')}{\partial n_{\vec r}} - G(\vec r,\vec r') \frac{\partial u(\vec r)}{\partial n_{\vec r}}\right) dS(\vec r).\end{multline}
Plugging in my new definition of $G(\vec r, \vec r')$, we see that any $u(\vec r)$ satisfying your PDE must satisfy
we can immediately see that any solution to the PDE must satisfy
$$ u(\vec r') = - \int_\Gamma G(\vec r,\vec r')g(\vec r) \ dS(\vec r) + c, \ \ \ \ \ \ (\ast )$$
where the constant $c$ is equal to $\frac 1 A \int_\Gamma u(\vec r) dS(\vec r)$, the average value of $u$ on $\Gamma$.

Having redefined the Green's function, I'll give you an explicit expression in the case where $\Omega$ is a two-dimensional circular disk of radius $1$. Here it is:
$$ G(\vec r, \vec r') = \tfrac 1 {2\pi} \ln | \vec  r -\vec r' | + \tfrac 1 {2\pi} \ln |\vec r - \vec r''|,$$
where $\vec r'' = \vec r'/|\vec r'|^2$ is the image of $r'$ under an inversion about the unit circle. It is similar to the Dirichlet Green's function, expect that we have a plus sign in front of the "image" term instead of a minus sign.
I believe that if you plug this Green's function into $(\ast )$ (with the constant $c$ chosen arbitrarily), you do indeed get a solution to your PDE, provided that $g$ obeys the consistency condition $\int_\Gamma g \ dS = 0$. This is stated in these lecture notes from my university, and also in Riley, Hobson and Bence, though I would love to see a rigorous proof! I wonder if anyone knows a good reference?
