As I wrote in the comments, I suspect that this question might have been prepared incorrectly: I don't see a method better than exploiting symmetry and integrating directly, and the resulting integral doesn't appear to have a closed form in terms of elementary functions, though it does have an unpleasant expression in terms of dilogarithms.
Since the integrand $\nabla f \cdot {\bf n}$ is invariant under symmetries of the icosahedron $S$, we need only integrate over a suitable triangle $T$ contained in a face $F$:
$$\iint_S \nabla f \cdot {\bf n} \,d\sigma = 20 \iint_F \nabla f \cdot {\bf n} \,d\sigma = 120 \iint_T \nabla f \cdot {\bf n} \,d\sigma .$$
Now, using the fact that the central angle between two adjacent vertices of $S$ is $\arccos \frac{1}{\sqrt 5}$ gives that the distance between the center $U$ of $F$ and a corner of $F$ is $\lambda R$, where $R$ is the radius of the sphere inscribed in $S$, and $\lambda := 3 - \sqrt{5} \approx 0.76393\!\ldots$.
In spherical coordinates centered at the center of $S$, $\nabla f = \rho^{-1} {\bf \hat \rho}$, so at a point $x$ on $F$ at a distance $r$ from $U$, $$(\nabla f)(x) = \frac{1}{\sqrt{R^2 + r^2}} \hat\rho ,$$ and then some straightforward trigonometery shows that $$\nabla f \cdot {\bf n} = \frac{R}{R^2 + r^2} .$$
Thus, in polar coordinates $(r, \theta)$ on $F$ centered at $U$, we have
$$\iint_E \nabla f \cdot {\bf n} \,d\sigma
= \int_0^{\frac{\pi}{3}} \int_0^{\frac{1}{2} R \lambda \sec \theta} \underbrace{\frac{R}{R^2 + r^2}}_{\nabla f \cdot {\bf n}} \underbrace{r \,dr \,d\theta}_{d\sigma}.$$
Computing the inner integral gives
$$\frac{1}{2} R \int_0^{\frac{\pi}{3}} \log\left(\frac{1}{4} \lambda^2 \sec^2 \theta + 1\right) \,d\theta .$$
This integral, however, does not appear to have a closed form in terms of elementary functions. It does have an explicit formula in terms of dilogarithms, but it's too complicated to be worth reproducing here. I would anyway be pleased to see a more tractable expression for this quantity.
Edit I made some progress on this integral---replacing $\frac{1}{2} \lambda$ with $\alpha$ and applying Feynman's trick we can reduce computing this integral to understanding $\frac{\arctan \beta \,d\beta}{\beta^2 - 3}$. Maple gives a slightly nicer expression, but it still involves dilogarithms. I didn't push the (rather ugly) computation to its end, but the form of the expression suggests that the integral can be written in terms of elementary functions and three instances of Legendre's chi function. It's not clear whether one can improve on this.
At any rate, the flux of $\nabla f$ across $S$ is thus
$$\iint_S \nabla f \cdot {\bf n} \,d\sigma = 60 \left[\int_0^{\frac{\pi}{3}} \log\left(\frac{1}{4} \lambda^2 \sec^2 \theta + 1\right) \,d\theta\right] R .$$ Integrating numerically gives that the quantity in brackets (just as well, the flux across the icosahedron with $R = 1$) is $13.35031\!\ldots$.
Since $\nabla \cdot \nabla f > 0$ everywhere, the Divergence Theorem gives that the flux across $S$ is between the fluxes of the spheres inscribing and circumscribing $S$, giving the cheap bounds
$$4 \pi R < \iint_S \nabla f \cdot {\bf n} \,d\sigma < 4 \pi \sqrt{R^2 + \lambda^2} .$$
Remark Another reason to suspect that the question is not what the preparer intended is that the analogous problem in dimension $2$---computing the flux of the vector field $\nabla \log\sqrt{x^2 + y^2}$ across a regular polygon in $\Bbb R^2$ with the origin on its interior---does admit an efficient solution that uses the Divergence Theorem and the fact that $\nabla \cdot \nabla \log\sqrt{x^2 + y^2} = 0$. A vector field in $\Bbb R^3 \setminus \{0\}$ with this behavior is $g(x, y, z) = (x^2 + y^2 + z^2)^{-1 / 2}$.
Best Answer
A parametrization of the disk S is $$S:\; (u,v) \rightarrow \; \left(0, v\cos(u), v\sin(u)\right) \;\;\text{with}\;\; 0\leq u\leq 2\pi, \; 0\leq v \leq \sqrt\frac{A}{\pi}.$$ Then $$\textbf{S}_{u} \times \textbf{S}_{v}=\left(0, -v\sin(u), v\cos(u)\right) \times \left(0, \cos(u), \sin(u)\right)=(-v,0,0)\implies \textbf{n}=(-1,0,0).$$
Actually, you don't need it, because $$\textbf{J} = -\sigma\nabla(V(1-x/L))=-\sigma(-V/L,0,0)$$ is constant and therefore $$\iint_S \textbf{J}\cdot\hat{\textbf{n}}\; dS=-\sigma(-V/L,0,0)\cdot(-1,0,0)A=\frac{\sigma VA}{L}.$$