You can see fairly quickly that an exact answer to this question is not going to be feasible because your functional transformation is to take the square root of $\sigma_t^2$, and the square root function has a countably infinite number of derivatives. This implies that a Taylor expansion is going to leave us with a countably infinite number of terms, most of which won't vanish.
So is all hope lost? No. Interestingly enough, for an appropriate choice of constant to expand around, the first order term in the Taylor expansion drops out, and the second term is known, so we can get a fairly good approximation.
Specifically, define $c_0 = \mathbb{E} \sigma_t^2$, which, as you point out in the question, is a known quantity. Note also that $\mathbb{V} \sigma_t^2$ is a known quantity (it is a reasonably complicated function of the parameters of the GARCH model - I can't remember the reference off the top of my head but you should be able to google it fairly easily). Using Taylor's theorem, we have:
\begin{equation}
\sigma_t = \sqrt{c_0} + \frac{1}{2} c_0^{-\frac{1}{2}}(\sigma_t^2 - c_0) - \frac{1}{8} c_0^{\frac{-3}{2}}(\sigma_t^2 - c_0)^2 + R
\end{equation}
where $R$ is the remainder term from the expansion. Taking expectations of both sides, you can see immediately that the first order term is going to vanish, since we will get $(\mathbb{E} \sigma_t^2 - c_0)$ which by definition of $c_0$ will equal zero. Further, note the bracketed portion of the second order term will become $\mathbb{E} (\sigma_t^2 - c_0)^2$ which again, given our initial choice of $c_0$, will simplify to $\mathbb{V} \sigma_t^2$. Thus, assuming that $\mathbb{E} R$ is small, we have:
\begin{equation}
\mathbb{E} \sigma_t \approx \sqrt{\mathbb{E} \sigma_t^2} - \frac{1}{8} (\mathbb{E} \sigma_t^2)^{\frac{-3}{2}} (\mathbb{V} \sigma_t^2)
\end{equation}
As I said before, $\mathbb{E} \sigma_t^2$ and $\mathbb{V} \sigma_t^2$ are both known expressions w.r.t. to the parameters of the model, and so the right hand side of this equation is known (albeit it would be very messy if I wrote it out with all the parameters from a GARCH model).
So how good is the above approximation? Initially I was tempted to suggest that it might be quite good (I was short on time and it seemed an easy way to tie off the answer). However, the comment made by @ZacharyBlumenfeld is a good one, and now I'm not so sure.
There are many factors at play here, some promising, and some not.
1) The series has alternating signs, i.e. $x_1 + x_2 - x_3 + x_4 - x_5 + ...$. This is good, as the conditions for a series with alternating signs to converge are much weaker than for a series with the same sign.
2) The Taylor coefficients get small quite quickly. They go up in the factorial, i.e. $\frac{1}{2!}, \frac{1}{3!}, \frac{1}{4!}$.
3) If $\omega > 1 - \alpha - \beta$, then $c_0 > 1$ and so $c_0^{k}$ gets small quickly as $k \rightarrow -\infty$. This is good. However, note that if $\omega < 1 - \alpha - \beta$ we get the opposite effect, which would be bad. So the parameters of the model matter.
4) We have moments of $\sigma_t^2$ increasing as the terms increase. This is the one that worries me the most. I'm fairly sure I read a paper a while back on laws of large numbers for GARCH processes that shows that under quite reasonable conditions, higher moments of $\sigma_t^2$ may not exist, i.e. are infinite. This would be very bad for our approximation.
If this question is really important to you, I would recommend leveraging the wonderful amount of computing power we have these days and running some simulations for different input parameter values and seeing how good the above approximation really is. Feel free to report back here with your findings!
Cheers, hope this helps.
-colin
