Here's an inelegant but straightforward, elementary solution. It begins with the case $n=1$.
Generally, a non-central chi-squared distribution arises as the sum of squares of independent Normal distributions, each with unit variance, but with varying and possibly non-zero means. (When all means are zero, we obtain the usual chi-squared distribution.) For $n=1$ there's just one such Normal distribution involved, which is convenient to express as the sum of its mean $\mu$ and a standard Normal variable $X$. Thus, we are concerned with finding the variance of $Y=(X+\mu)^2$. By definition, this is
$$\operatorname{Var}(Y) = E[Y^2] - E[Y]^2 = E[(X+\mu)^4] - E[(X+\mu)^2]^2.$$
We can read these expectations directly off the moment generating function for $X+\mu$, which is
$$\psi_{X+\mu}(t) = e^{\mu t} \psi_X(t) = e^{\mu t} e^{t^2/2} = e^{\mu t + t^2/2}$$
as demonstrated at https://stats.stackexchange.com/a/176814/919, for instance. The Maclaurin series for $\psi$ is
$$\psi_X(t) = 1 + \mu t + \frac{1+\mu^2}{2!}t^2 + \frac{3\mu+\mu^3}{3!}t^3 + \frac{3+6\mu^2+\mu^4}{4!}t^4 + \cdots.$$
The numerators of the fractions are the corresponding moments, whence
$$E[(X+\mu)^4] = 3+6\mu^2+\mu^4$$
and
$$E[(X+\mu)^2]^2 = (1 + \mu^2)^2,$$
yielding
$$\operatorname{Var}(Y) = 3+6\mu^2+\mu^4 - (1 + \mu^2)^2 = 2(1 + 2\mu^2).\tag{1}$$
The hard part has been done, because the general case $n\ge 1$ reduces to this one.
The underlying geometric idea is that we can rotate the coordinate system so that the first coordinate is parallel to the vector of means $(\mu_1,\mu_2,\ldots,\mu_n)$ and the other coordinates are orthogonal to it. Because a rotation is a linear transformation, the coordinates continue to have a jointly multinormal distribution. Because rotations are orthogonal, the new variables continue to have unit variances. All the "noncentrality," however, has been located in the first coordinate alone because the remaining variables are standard Normal (their means are all zero).
These contours represent the underlying multivariate Normal distributions. The non-central chi-squared distribution is that of their squared distance to the origin. The distributions are the same in both panels, because a rotation around the origin does not change the distances to it. However, only the first coordinate in the rotated (right hand) panel has a nonzero mean; the remaining coordinates have standard Normal distributions.
I'll illustrate by doing the algebra for $n=2$. Let's index $Y$, $X$, and $\mu$ with $i=1,2$, where the $X_i$ are independent standard normal variables and $Y_i=\mu_i+X_i$. The variable $Y=Y_1^2 + Y_2^2$ has a noncentral chi-squared distribution with noncentrality parameter $\mu_1^2 + \mu_2^2$. Indeed, writing $\mu^2 = \mu_1^2 + \mu_2^2$ and assuming it is nonzero, we may perform a little algebra (completing the square) to obtain
$$\eqalign{
Y &= (\mu_1+X_1)^2 + (\mu_2+X_2)^2 \\
&= (\mu_1^2+\mu_2^2) + 2\mu_1X_1 + 2\mu_2X_2 + X_1^2 + X_2^2\\
&= \left(\mu + \frac{\mu_1}{\mu}X_1 + \frac{\mu_2}{\mu}X_2\right)^2 + \left(\frac{\mu_2}{\mu}X_1 - \frac{\mu_1}{\mu}X_2\right)^2 \\
&=(\mu + Z_1)^2 + Z_2^2
}$$
where
$$Z_1 = \frac{\mu_1}{\mu}X_1 + \frac{\mu_2}{\mu}X_2$$
and
$$Z_2 = \frac{\mu_2}{\mu}X_1 - \frac{\mu_1}{\mu}X_2.$$
The $Z_i$ are linear combinations of independent standard Normal variables. Their means therefore are $0$ and their variances are the sum of squares of the coefficients, both of which are
$$\operatorname{Var}(Z_i) = \left(\pm\frac{\mu_1}{\mu}\right)^2 + \left(\frac{\mu_2}{\mu}\right)^2 = \frac{\mu_1^2 + \mu_2^2}{\mu^2}=1.$$
Moreover, a similar calculation shows the $Z_i$ are uncorrelated. Since they are jointly Normal, they must be independent. Accordingly,
$$\operatorname{Var}(Y) = \operatorname{Var}((\mu + Z_1)^2) + \operatorname{Var}(Z_2^2).$$
Our previous formula $(1)$ applies directly to both terms (the noncentrality parameter in the second term is $0$), giving
$$\operatorname{Var}(Y) = 2(1+2\mu^2) + 2(1 + 2(0^2)) = 2(2 + 2\mu^2).$$
The formula for general $n$ works out the same way (use induction or apply linear algebra), giving the general formula
$$\operatorname{Var}(Y) = 2(n + 2\mu^2)$$
where now $\mu^2 = \mu_1^2 + \mu_2^2 + \cdots + \mu_n^2$ is the noncentrality parameter.
