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!