When we prove the average property of harmonic function, we use a formula
\begin{align} & \int_{B_r(x)}\triangle u\,dy=\int_{B_r(x)}\text{div}(\triangledown u)\,dy \\[6pt] = {} & \int_{\partial B_r(x)}\triangledown u\cdot v\,dS \\[6pt] = {} & \int_{\partial B_r(x)}\frac{\partial u}{\partial r}\,dS \tag{$*$} \\[6pt] = {} & r^{n-1}\frac{d}{dr}\int_{\partial B_1(0)}u(x+rw)\,dS \end{align}
I want to know what $\dfrac{\partial u}{\partial r}$ exactly mean, and why the equality $(*)$ is true. Thank you for your help.
The divergence theorem says that $$\int_{B_r(p)}\Delta u(x)\>{\rm d}(x)=\int_{B_r(p)}{\rm div}\bigl(\nabla u\bigr)(x) \>{\rm d}(x)=\int_{\partial B_r(p)}\nabla u(x)\cdot{\bf n}(x)\>{\rm d}\omega\ ,$$ where ${\bf n}(x)$ denotes the outward unit normal and ${\rm d}\omega$ the surface element on $\partial B_r(p)$.
Now $\nabla u(x)\cdot{\bf n}(x)$ is the directional derivative of $u$ in direction ${\bf n}$: $$\nabla u(x)\cdot{\bf n}(x)=\lim_{t\to0+}{f(x+t{\bf n})-f(x)\over t}\ ,$$ and it is customary to denote this "normal derivative" by ${\displaystyle{\partial u\over\partial n}}$. When the surface in question is a sphere centered at the origin (which, strictly speaking, is not the case here) we may write ${\partial u\over \partial r}$ instead of ${\partial u\over\partial n}$, because this directional derivative coincides with the partial derivative ${\partial u\over \partial r}$ when the cartesian coordinates are replaced by spherical coordinates.