I don't know of a fully intuitive derivation, but there are some informal arguments that give the circular law with a relatively small amount of calculation.

Let $M$ be a matrix where the entries are iid with mean zero and variance one. One can begin with the determinant formula

$$ \log |\det( M - z )| = \sum_{j=1}^n \log |\lambda_j - z|.$$

The circular law suggests that the eigenvalues $\lambda_j$ should be uniformly distributed in the disk of radius $\sqrt{n}$, so one should be proving something like

$$ \log |\det( M - z )| \approx \frac{1}{\pi} \int_{|w| \leq \sqrt{n}} \log |w-z|\ dw.$$

A routine calculation (e.g. using Jensen's formula, or the fundamental solution for the Laplacian) reveals that the RHS is equal to $n \log |z|$ when $|z| \geq \sqrt{n}$ and $\frac{1}{2} n \log n - \frac{1}{2} n + \frac{1}{2} |z|^2$ for $|z| \leq \sqrt{n}$. So heuristically, the circular law is equivalent to the approximations

$$ |\det(M-z)| \approx |z|^n$$

for $|z| \geq \sqrt{n}$ and

$$ |\det(M-z)| \approx n^{n/2} e^{-n/2} e^{|z|^2/2}$$

for $|z| \leq \sqrt{n}$. Here one should interpret the $\approx$ symbol rather loosely (in particular, polynomial factors in $n$ should be considered negligible).

However, by the Leibniz formula for determinants and the iid mean zero variance one nature of the entries (which makes all the covariances between the terms in the Leibniz formula vanish), one can easily compute that

$$ {\bf E} |\det(M-z)|^2 = \sum_{j=0}^n |z|^{2j} \frac{n!}{j!}.$$

(This type of calculation goes back to an old paper of Turan.)

For $|z| \gg \sqrt{n}$, the $|z|^{2n}$ term on the RHS dominates, while for $|z| \ll \sqrt{n}$, the RHS is most of the Taylor series for $n! e^{|z|^2}$. The claim then morally follows from Stirling's approximation.

(Incidentally, my recent paper with Van Vu on local versions of the circular law basically proceeds by making the above argument rigorous; see also recent work of Bourgarde, Yau, and Yin. The idea of controlling the spectrum of an iid matrix through its log-determinant goes back to the early work of Girko.)

--

Note also that without too much calculation, one can see that the limiting law of the spectrum should be invariant with respect to rotations around the origin (especially if one assumes that the entries of the iid matrix are similarly invariant, e.g. they are complex gaussian). From the matrix inequality $\sum_{j=1}^n |\lambda_j|^2 \leq \hbox{tr}(M M^*)$ and the law of large numbers we also see that the typical size of an eigenvalue $\lambda_j$ should be of the order of $\sqrt{n}$. These facts fall well short of the full circular law but are certainly consistent with that law.