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-1} \nabla u(x)\cdot h\, dx =\left(\int_\Omega p|\nabla u(x)\cdot e|^{p-1} \nabla u(x)\,dx\right)\cdot h $$$$ 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-1} \nabla u(x)\,dx $$$$ 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-1} \nabla u(x)\,dx\right|}$$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-1}\nabla u$$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=1$$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!