Ok I had a question I think I can almost answer it but I miss one step:
Let $\partial M$ be a closed surface in $\mathbb{R}^3$, $x_0 \in \partial M$ than show this limit: $$\lim_{\substack{x\rightarrow x_0 \\ x\in M}} \int_{\partial M} \frac{1}{||y-x||} n_y \cdot \nabla_y \frac{-1}{||y-x||} dS_y = \lim_{\substack{x\rightarrow x_0 \\ x\in M}} \int_{\partial M} \frac{n_y \cdot (y-x)}{||y-x||^4} dS_y = \infty$$
The condition on $\partial M$ is free to choose. So suppose that it is smooth as much as you need. But ideally I would like to prove it for piece wise $C^1$ surface.
So has anybody idea how to show this limit?
Edit: Ok my best shot when $\partial M$ is convex. But I can't prove it yet, I'm experiencing the same problem as when I show that $\int_0^1 \frac{1}{t^\alpha} dt$ diverge(for $\alpha > 1 $) by this technique:
$$ \int_0^1 \frac{1}{t^\alpha} dt \geq \int_0^x \frac{1}{t^\alpha} dt \geq \frac{x}{x^\alpha} \rightarrow \infty$$
as $x$ goes to zero. But IT doesn't work for $\alpha =1$ yet the integral diverge.
So here goes my try, I need two lemmas about convex surfaces.
Lemma1: It is bound on how big can be angle between $n_y$ and $y-x$.
Let $\partial M$ is convex surface in $\mathbb{R}^n$. Than for every $x$ in interior of $\partial M$ ie $x\in M$. This inequality holds: $$ \forall y\in \partial M: n_y\cdot \frac{y-x}{\|y-x\|} \geq \frac{\min_{z\in \partial M} \|x-z\|}{\max_{z\in \partial M} \|x-z\|}\geq \frac{\text{dist}(x,\partial M)}{\text{diam}(\partial M)}$$ $n_y$ is outer normal at point $y$
This lemma is quite straight forward just draw a picture in 2d and you'll see it.
Lemma2: This is bound on how small can get and surface area of convex surface at vicinity of some point $x$
Let $\partial M$ is convex surface in $\mathbb{R}^n$ and $x$ is point in interior of $\partial M$ ie $x\in M$. Denote $\rho = \text{dist}(x,\partial M) = \text{dist}(x,x_0)$ for some $x_0 \in \partial M$. Than: $$\int_{\partial M \cap B(x_0,2\rho)} 1 \, dS \geq \alpha_{n-1} \left(\frac{\sqrt{3}}{2} \rho \right)^{n-1} $$ where $\alpha_n$ is volume of unit ball in $\mathbb{R}^n$.
Again draw this in 2d and you'll see it. Maybe I will add pictures if I find out how.
Now to the limit:"
denote $\rho = \text{dist}(x,\partial M)=\|x-x_0\|$, $D = \text{diam}(\partial M)$
$$ \int_{\partial M} \frac{1}{||y-x||} n_y \cdot \nabla_y \frac{-1}{||y-x||} dS_y = \int_{\partial M} \frac{n_y \cdot (y-x)}{||y-x||^4} dS_y \geq$$
$$ \int_{\partial M} \frac{n_y \cdot (y-x)}{||y-x||^4} dS_y \geq \int_{\partial M} \frac{\rho \|y-x\|}{D\|y-x\|^4} dS_y \geq $$
$$ \int_{\partial M \cap B(x_0,2\rho)} \frac{\rho \|y-x\|}{D\|y-x\|^4} dS_y \geq \int_{\partial M \cap B(x_0,2\rho)} \frac{\rho^2}{D 2^4 \rho^4} dS_y \geq$$
$$ \frac{\rho^2}{D 2^4 \rho^4} \int_{\partial M \cap B(x_0,2\rho)} dS_y \geq \frac{\rho^2}{D 2^4 \rho^4} \pi \left(\frac{\sqrt{3}}{2} \rho \right)^{2} = \frac{1}{D 2^4 } \pi \frac{3}{4}$$
But that does not got to infinity :((
I have one minor interesting observation. Let $\Gamma$ be any surface in $\mathbb{R}^3$. $\Gamma$ can be self intersecting for example: embedded Klein bottle in $\mathbb{R}^3$.
Than I think that
$$-\frac{1}{4\pi}\int_{\partial M} n_y \cdot \nabla_y \frac{1}{||y-x||} dS_y$$
is something like winding number of $\Gamma$ around $x$. Can anybody elaborate on that? What number would I get if $\Gamma$ would be Klein bottle? I think that Klein bottle doesn't have something like interior so is the integral zero everywhere? Therefore is it possible in 2D to have closed, regular $C^1$ curve that has empty interior?
Denote $\Gamma(x)=\frac1{|x|}$. For fixed $x\in M$ as a function of $y$ it is harmonic in $\mathbb R^3\backslash \bar M$. Integrating by parts yields $$ \int_{\mathbb R^3\backslash M}|\nabla_y\Gamma(x-y)|^2\,dy= \int_{\mathbb R^3\backslash M}\nabla_y\Gamma(x-y)\cdot\nabla_y\Gamma(x-y)\,dy= $$ $$ -\int_{\mathbb R^3\backslash M}\Gamma(x-y)\Delta \Gamma(x-y)\,dy+ \int_{\partial M}\Gamma(x-y)n_y\cdot\nabla_y\Gamma(x-y)\,dy= $$ $$ =\int_{\partial M}\Gamma(x-y)n_y\cdot\nabla_y\Gamma(x-y)\,dy $$ since $\Gamma\nabla_y\Gamma$ is vanishing at infinity fast enough, as $|y|^{-3}$. Here $n_y$ is an inner unit normal to $M$. So the integral in question is equal to (may be up to the sign, depending upon choice of $n_y$) $$ I(x)=\int_{\mathbb R^3\backslash M}|\nabla_y\Gamma(x-y)|^2\,dy. $$
Denote $d=d(x)=\mathrm{dist}(\partial M,x)$, $r=|x-y|$ and $B_d(x)$ a ball of radius $d$ centered at $x$. To get a grasp on the last integral near the boundary let's take an inversion relatively the ball's surface $S_d(x)$, $r\to \frac{d^2}r$. Let $\tilde M_x$ be the image of ${\mathbb R^3\backslash M}$ under the inversion.
Substituting this transform in the integral gives $$ I(x)=\int_{\mathbb R^3\backslash M}|\nabla_y\Gamma(x-y)|^2\,dy= \int_{\mathbb R^3\backslash M}r^{-4}\,dy= $$ $$ \int_{\mathbb R^3\backslash M}r^{-4} r^2\,\cos\theta\,d\theta d \varphi dr= $$ $$ \int_{\tilde M_x}\frac{1}{\frac{d^4}{r^2}}\cos\theta\, d\theta d \varphi\, \frac{d^2}{r^2}dr= $$ $$ \frac1{d^2}\int_{\tilde M_x}\cos\theta\, d\theta d \varphi dr= \frac1{d^2}\int_{\tilde M_x}\frac1{r^2}\,dy. $$ As is seen from the formula in the spherical coordinates the integral sums up solid angles $A(r)$ of intersections $\tilde M_x\cap\{r=\mathrm{const}\}$ as seen from $x$: $$ I(x)=\frac1{d^2}\int_0^d A(r)\,dr. $$
Roughly speaking, the more space $\tilde M_x$ is taking in the ball $B_d(x)$ the larger is $I(x)$. So for $I(x)$ to be blowing up as $d\to0$ $\tilde M_x$ should take large enough part of $B_d(x)$. If $x$ is near the boundary $\partial M$ then the largest contribution to $I(x)$ comes from the part of $\mathbb R^3\backslash M$ that is near to $x$.
Let's consider the case $\partial M\in C^1$. Then for small $d$ sphere $S_d(x)$ touching $\partial M$ looks in vicinity of $x$ as if it is touching a plane. The compliment of $M$ is looking like a half-space.
An inversion of a tangent half-space to a unit sphere is a ball of radius $1/2$ touching the sphere from inside. It is easy to see for the planar case as the tangent line turns into a circle and in the three-dimensional case there is a rotational symmetry. The integral for $I(x)$ can be evaluated explicitly here, since the equation of the ball in the spherical coordinates is $0<r<\sin\theta$, $0<\theta<\pi/2$, and $$ A(r)={\int_0^{2 \pi } \left(\int_{\arcsin r}^{\pi/2} \cos \theta \, d\theta \right) \, d\varphi }=2\pi(1-r), $$ $$ \int_0^1 A(r)\,dr=\pi. $$ For a ball of radius $d$ the answer is multiplied by $d$, so $$ \frac1{d^2}\int_0^d A(r)\,dr=\frac1{d^2}\pi d=\frac\pi{d}. $$ Hence for bounded domains from $C^1$ the asymptotic is $I(x)\sim\frac{\pi} d$ as $d\to0^+$. In particular $\lim_{x\to x_0}I(x)=\infty$.
For a convex domain the answer is evidently the same and for polyhedrons it should be too. For the class of piecewise $C^1$ domains I think the answer is negative if $x_0$ is a point of an inward spike (or a cusp). In such cases $A(r)$ could be small enough to get $\int_0^d A(r)\,dr\le Cd^2$ for some points $x_n$ s.t. $x_n\to x_0$.
Note. The condition for $\lim_{x\to x_0}I(x)=\infty$ written above seems alike conditions for the boundary point $x_0$ to be regular. Naturally the question arises if it exactly so. If yes, it would give some known sufficient conditions for $x_0$ to have the required property, such as the exterior sphere condition or the exterior cone condition.