Let $\Omega\subset \Bbb R^d$ be a bounded $C^1$ domain. Let $u:\Bbb R^d\to \Bbb R$ be a function in $C^2_b(\Bbb R^d)$. I would like to compute the following limit: for $x\in \partial \Omega$
$$L= \lim_{s\to 1}(1-s)\int_{\Omega}\frac{(u(x)-u(y))}{|x-y|^{d+2s}} d y. $$
Here is what I did so far:
Let $r>0$ be arbitrarily small enough. Then as $u$ is bounded, we have
$$\begin{align}&\lim_{s\to 1}(1-s)\int_{\Omega\cap \{|x-y|\geq r\}}\frac{|u(x)-u(y)|}{|x-y|^{d+2s}}dy\\&\leq C \lim_{s\to 1}(1-s)\int_{|x-y|\geq r}\frac{dy}{|x-y|^{d+2s}} \\&= Cc_d\lim_{s\to 1}(1-s) \int_r^\infty t^{-2s-1} dt= 0. \end{align}$$
so that
$$L= \lim_{s\to 1}(1-s)\int_{\Omega \cap B_r(x)}\frac{(u(x)-u(y))}{|x-y|^{d+2s}} d y. $$
Using the fundamental theorem of calculus and the relation $\nabla [|x|^\alpha]= \alpha x|x|^{\alpha-2}$,
$$\begin{align}&(1-s)\int_{\Omega \cap B_r(x)}\frac{(u(x)-u(y))}{|x-y|^{d+2s}} d y\\ & = -\int_0^1 dt (1-s)\int_{\Omega \cap B_r(x)}\frac{\nabla u(x+ t(y-x))\cdot (y-x)}{|x-y|^{d+2s}}dy\\&= \int_0^1 dt\frac{ (1-s)}{d-2(1-s))}\int_{\Omega \cap B_r(x)} \nabla u(x+ t(y-x))\cdot \nabla_y [|x-y|^{-d+2(1-s)}]dy \end{align}$$
Therefore my initial question could be resumed to computaion of
$$\lim_{s\to 1} (1-s)\int_{\Omega \cap B_r(x)} \nabla u(x+ t(y-x))\cdot \nabla_y [|x-y|^{-d+2(1-s)}]dy.$$
Any idea on to move further?
My feeling is the should be a multiple factor of $\frac{\partial u}{\partial n}(x)= \nabla u(x).n(x)$ where $n(x)$ is the normal derivative on $\partial \Omega$ at the point $x$. Don't be mind corrupted, I may have wrong expectation.
Here's a sketch. I have not finished the last step yet, but maybe this will be useful anyhow... I will replace $d$ in your notation by $n$ to avoid confusion with derivatives and differentials. Also, the exponent of the denominator should be $n+s$ for this question to make sense, and the limit must be one sided from below.
The key idea is that since the quantity of interest $L$ is local we can transfer the problem to a space which is locally the same, but easier to understand (a subset of the half space). Since the domain is $C^1$, there is a $C^1$ chart $\phi$ mapping a neighborhood $U\subset H_n$ of the origin diffeomorphically to a neighborhood of $x$, where $H_n$ is an n dimensional half space, and such that $\phi(0)=x$.
By your first observation we can take the integral to be over a small ball $B_r$ of radius $r$ near $x$, small enough that it lies in the image $\phi(U)$. Then we can do a change of variables to rewrite $L$ as $$L=\lim_{s\rightarrow 1} (1-s)\int_{\phi^{-1}(B_r)} \frac{u(\phi(0))-u(\phi(v))}{|\phi(0)-\phi(v)|^{n+s}} J(v) dv$$ where $J(v)=|D\phi(v)|$ is the Jacobian of $\phi$. We can restrict the domain of this new integral in the same way as before, because under $\phi$ the exterior of any ball around the origin is mapped to the exterior of a ball about $x$ in the original space (and the integral over such a region vanishes in the limit). So replace the domain of integration by a new ball $B_R\cap H_n\subset \phi^{-1}(B_r)$.
We approximate $$u(\phi(0))-u(\phi(v)) = -D(u\circ \phi)(0) V + \mathcal{o}(v) = -Du(x)D\phi(0)v + \mathcal{o}(v)$$ $$J(v) = J(0) + \mathcal{o}(v)$$ and lastly $$\phi(0)-\phi(v) = -D\phi(0)v + \eta(v)$$ with $\eta(v) = \mathcal{o}(v)$. We give a name to this last error term, because it will require some extra consideration.
Writing $v=t w$ with $t=|v|$ and $w\in S^{n-1}\cap H_n$, we transform the integral to spherical coordinates:
$$ L= \lim_{s\rightarrow 1}(1-s) \int_{S^{n-1}\cap H_n} dw \int_0^R dt \left[ \frac{-Du(x) D\phi(0)w t^n J(0) +\mathcal{o}(t) t^{n-1} }{|D\phi(0) w t + \eta(wt)|^{n+s}} \right]$$
We want to show that the $\mathcal{o}(v)$ terms can be ignored, and evaluate what remains. We estimate the integral in several steps: