Why isn't the Welch–Satterthwaite equation nonsensical? As I understand it, the Welch–Satterthwaite equation says that if $s_i^2$ is the sample variance of group $i$ (with $n_i$ samples), then, assuming that the measurements in each group are iid's that are normally distributed, the statistic $\sum a_i s_i^2$ has roughly $\chi^2$ distribution with degrees of freedom $\frac{(\sum a_i s_i^2)^2}{\sum \frac{(a_is_i^2)^2}{n_i-1}}$.
This equation doesn't take into account the $\sigma_i$'s, and that makes me very suspicious...!
Here's a quick example to show how weird this formula is. Assume for a second that we apply this equation to the case that $a_1=1$ and $a_i=0$ for $i>1$. Then this equation says that $s_1^2$ has roughly $\chi^2$ distribution with $n_1-1$ degrees of freedom. But that's incredibly false! What if $\sigma_1^2$ is huge?
When you take more than one $a_i$ to be nonzero, more ridiculousness ensues.
So what, if anything, am I getting wrong?
 A: It often happens, in practical applications, that two or more sample variances are available and it is desired to add them.  For simplicity, assume just two sample variances and assume they are as you describe in your first sentence.  Also assume the respective population variances are unknown.  (If they were known, there would be no reason to bother with the sample variances.) Each sample variance is gamma distributed:
 $$s_i^2 \thicksim \Gamma(k_i,\theta_i)$$
where $k_i$ are the respective shape parameters, equal to $\nu_i$/2, and $\theta_i$ = 2$\sigma_i^2/\nu_i$ are the respective scale parameters.
The sum will be gamma distributed iff the two scale parameters are equal.  Since the population parameters are unknown, so are the scale parameters.  Hence the sum is either a sample from a finite mixture of gamma distributions, as @whuber elegantly showed or a sample from an infinite sum of gamma distributions, as per Moschopoulos’s paper (P.G. Moschopoulos, The Distribution of the Sum of Independent Gamma Random Variates, Ann. Inst. Statist. Math. 37 (1985), Part A, 541-544). The latter would likely be the case.  These are very interesting, but not so useful in practical applications.
So the Welch-Satterthwaite approximation consists of defining the sum of the two sample variances as being a sample from a fictitious gamma distribution, with its shape and scale parameters computed from those of the two constituent sample variance gamma distributions. But you do not know the two population variances, so you substitute the respective sample variances. Then degrees of freedom is computed from one of several equations: I have seen at least 3 inequivalent equations for this.
So does it work?  Yes and no. The intention is to obtain a gamma distributed sum, so the square root would be $\chi$ distributed, like our customary sample standard deviations, and then it would be feasible to use critical t values, construct confidence intervals, etc. The approximation fails if the sum is essentially just the larger summand. I have seen a paper where one population variance is known, and used in the Welch-Satterthwaite approximation, along with the sample variance for the other variate. I find this to be puzzling and the associated degrees of freedom can be rather large. However, used reasonably, with sample variances that are not far different, the Welch-Satterthwaite approximation has some utility.
For more information, go here: http://www-personal.umd.umich.edu/~fmassey/gammaRV/
and click on subsection 4.1 under the Welch-Satterthwaite approximation.  This will download a 3 page document that shows how the approximation arises.
A: 
This equation doesn't take into account the $\sigma_i$s, and that makes me very suspicious...!

You are right to be suspicious, because the formulation you have expressed is wrong.  The Welch-Satterwaite approximation should include a scaling term to give the chi-squared distribution, so it should be:
$$\frac{\sum_i a_i s_i^2}{\sum_i a_i \sigma_i^2} 
\overset{\text{Approx}}{\sim} \frac{\chi^2_{DF}}{DF} 
\quad \quad \quad \quad \quad
DF = \frac{(\sum_i a_i s_i^2)^2}{\sum_i \frac{(a_is_i^2)^2}{n_i-1}}.$$
The version you have written in your question lacks the scaling term $\sum_i a_i \sigma_i^2$ on the left-hand-side to scale the expected value to a unit value, and the scaling term $DF$ on the right-hand-side to give a corresponding scaled chi-square (with unit expectation).  (Obviously you can shift the scaling term on the right-hand-side to the left via multiplication if you prefer.)  Using this correct formulation of the approximation, in the special case where you have $a_1=1$ and all other $a_i=0$ you get the result:
$$\frac{s_1^2}{\sigma_1^2} 
\overset{\text{Approx}}{\sim} \frac{\chi^2_{n_1-1}}{n_1-1} .$$
This approximation holds exactly in the case of normal data, and is a reasonable large-sample approximation for non-normal data with finite kurtosis.  (In the cases where the data is not mesokurtic you can adjust the degrees-of-freedom using an alternative formula in O'Neill 2014.)
