Sums of the form $S_k=\sum_{n\geq 1}\frac{1}{\sinh^{2k}(n \pi)}$ and the residue theorem. Let us define the sums

$$S_k=\sum_{n\geq 1}\frac{1}{\sinh^{2k}(n \pi)}$$

A few days ago, i was able to give a surprisingly simple derivation of the fact
that
$$
S_1=\sum_{n\geq 1}\frac{1}{\sinh^2(\pi n)}=\frac{1}{6}-\frac{1}{2\pi}
$$
as an answer to this nice old question in terms of contour integration.
Being curious i immediatly tried to generalize this result to higher orders of $\sinh$, which as it turned out will not work that easily. For example, considering the case of $k=2$ it turns out that  
$$\text{Res}\left(\frac{
\cot(\pi z)}{ \sinh^4(\pi z)},z=z_k\right)=
\begin{cases}
\frac{1}{\pi\sinh^4(\pi k)} \quad\text {if}\,\,k \in \mathbb{Z}/0,\\
-\frac{1}{\pi\sinh^4(\pi k)} -\frac{4}{3\pi\sinh^2(\pi k)} \quad\text {if}\,\,ik \in \mathbb{Z}/0
\end{cases}
$$
So integrating over a large quadratic contour gives us absolutly no information about $S_2$, since the residues cancel (For details of the exact procedure please have a look at the post linked above). This is in fact a simple consequence of the transformation properties of $\sinh^{2k}(z)$ along the imaginary axis: $\sinh^{2k}(i y)=(-1)^k\sin^{2k}(y)$ 
Looking at $S_3$ with the same approach, it turns out that it contains $S_2$ as well as $S_1$ so it is clear that we need this value to make any progress in a generalization at all. 
I also want to mention that a Laplace/Mellintransform approach seem to suffer from the same cancellations then what i tried above, so this  is maybe also not the way to go...:/
My questions are

1.) Is there any chance to proof the result for $S_2$ (and $S_k$ in general if possible) using an an approach which is close to the one i used to calculate $S_1$ ? 
2.) What is a general approach to derive calculate such sums using a minimum of knowledge about the realm of elliptic integrals, Jacobi theta functions etc. (where i'm sure they can be derived but i have far too less knowledge in this field to do this by my own) 

PS:
It turns out that a closed form is, according to mathematica, given by
$$
S_2=-\frac{11}{90}+\frac{1}{3\pi}+\frac{\bar{\omega}^4}{7680
\pi^4}
$$
where $\bar{\omega}$ is the Lemniscate constant.
 A: ($1$) IMHO there is little chance to compute $S_2$ using contour integration due to existence of the lemniscate constant in all $S_n, n\ge 2$ (see below).
($2$) Here is a general approach. The normalized Eisenstein series can be written as
$$G_{2k}(i)=\sum_{m,n\in\mathbb R, mn\not=0}\ \frac{1}{(m+ni)^{2k}}=2\zeta(2k)+2\sum_{n>0}\sum_{m\in\mathbb R}\frac{1}{(m+ni)^{2k}}$$
Where the inner sum can be computed through differentiating Mittag-Leffler expansion of $\cot$:
$$\sum_{m\in\mathbb R}\frac{1}{(m+ni)^{2k}}=\sum_{j=1}^{k}a_{k,j}\text{csch}^{2j}(n\pi)$$
Thus the following recurrence
$$G_{2k}(i)=2\zeta(2k)+2\sum_{j=1}^{k}a_{k,j}S_k$$
with well-known initial value $S_1=\frac{1}{6}-\frac{1}{2\pi}$ (see for instance J. Borwein's monograph Pi and AGM) allows us to compute all $S_n$ in terms of $G_{2k} (k=2,\cdots,n)$, all of which are expressible via products of $G_4(i)=\frac{\Gamma \left(\frac{1}{4}\right)^8}{960 \pi ^2}, G_6(i)=0$ (due to theory of modular forms). Therefore we conclude that $S_n$ lie in the algebra of $\mathbb Q$ generated by $\frac{1}{\pi}, \Gamma\left(\frac14\right)$, for instance
$$S_5=\frac{14797}{187110}-\frac{64}{315 \pi }-\frac{\Gamma \left(\frac{1}{4}\right)^{16}}{2580480 \pi ^{12}}-\frac{41 \Gamma \left(\frac{1}{4}\right)^8}{90720 \pi ^6}$$
