Optimizing the gradient norm on the unit sphere

I do not believe there is a simple formula to express $I(u)$, but for sure for most of the functions the inequality $$ I(u)<\int_\Omega |\nabla u|^p\, dx $$ is sharp. For example if $\Omega=B$ is a ball and $u(x)=f(|x|)$ is a radial function, then $\nabla u$ is a vector field orthogonal to the sphere $\mathbb{S}^{d-1}$. Then $|\nabla u(x)|=|\nabla u(x)\cdot e|$ only for points $x$ on the line passing through $0$ and parallel to $e$ for all other points the inequality is sharp.

To have equality you would have to have supremum under the sign of the integral, but that would be a different expression.

More of a comment than an answer, really, but too long for the comment box: For a fixed smooth function $u$ the map $$ e\mapsto U(e):=\int_\Omega |\nabla u(x)\cdot e|^p dx $$ is differentiable as a function of $e\in\mathbb R^{d}$, and its differential in the direction $h$ is simply $$ DU(e).h =\int_\Omega p|\nabla u(x)\cdot e|^{p-2} \nabla u(x)\cdot h\, dx =\left(\int_\Omega p|\nabla u(x)\cdot e|^{p-2} \nabla u(x)\,dx\right)\cdot h $$ A maximizer $e$ on the unit sphere must then satisfy the first-order optimality condition $ DU(e)\cdot h=0$ for all tangent directions $h\in T_e\mathbb S^{d-1}\Leftrightarrow h\cdot e=0$, which means here that $DU(e)$ must be colinear to $e$. In other words, $$ e=C \int_\Omega p|\nabla u(x)\cdot e|^{p-2} \nabla u(x)\,dx $$ for some normalization constant $C=\frac{\pm1}{\left|\int_\Omega p|\nabla u(x)\cdot e|^{p-2} \nabla u(x)\,dx\right|}$ (provided the denominator does not vanish). Of course we have a $\pm$ degree of freedom due to the invariance $U(e)=U(-e)$. I don't know how much one can exatract from this integral condition, but at least it is clear that the reasonable guess $e=C\int |\nabla u|^{p-2}\nabla u$ is too naive and does not work (since it does not satisfy a priori this integral condition).

Note that for $p=2$ the solution is obviously given by $e= C\int \nabla u$, the average gradient (provided it is not zero, of course), so the functional is somehow the "directional $TV$ norm" $I(u)=\int |\partial_e u|$ in the average (most varying) direction $e=C\int \nabla u$.

Interesting functional!

A bit of a longer comment and expansion on Piotr's answer: if you restrict to the case $p$ is even, you can write your integral as $$ \int |\nabla u \cdot e|^p = \int (\nabla u \otimes \cdots \otimes \nabla u) \cdot (e \otimes \cdots \otimes e) $$ where the tensor product is taken over $p$ times. Since the arguments are manifestly symmetric, you have $$ I(u) \leq \sup_{T \in \mathcal{S}^P\mathbb{R}^d, |T| = 1} \int (\otimes^p \nabla u) \cdot T =: J(u) $$ Here $\mathcal{S}^p\mathbb{R}^d$ is the set of symmetric rank $p$ tensors over $\mathbb{R}^d$, with induced inner product. The optimization of $J(u)$ is a linear problem and can be easily solved with $$ T = \frac{\int \otimes^p \nabla u}{| \int \otimes^p \nabla u|}$$ and $J(u) = | \int \otimes^p \nabla u |$.

If you are lucky enough that this $T$ is a pure tensor (of the form $e \otimes e \otimes \cdots \otimes e$) then $I(u) = J(u) = | \int \otimes^p \nabla u |$. However, for generic $u$ this cannot be expected, and you get $I(u) < J(u)$ in this case.

As a final note, you have that $$ \left| \int \otimes^p \nabla u \right| \leq \int |\otimes^p \nabla u| = \int |\nabla u|^p $$ The inequality step is strict as long as $\nabla u$ is not almost everywhere parallel to a fixed vector. So this method, while it doesn't give you the exact value of $I(u)$, can at times give you a slightly improved upper bound.