Let ${\bf X}$ be a $(2n,k)$ multinomial random variable, corresponding to the experiment of placing $2n$ balls inside $k$ urns, with uniform probability. Then
$$P_{\bf X} = \frac{(2n)!}{x_1!x_2! \cdots x_k!} (1/k)^{2n} \left[\sum x_i = 2n\right] \tag 1$$
Let $E$ be the event that all urns have an even number of balls. Then
$$k^{2n} \, P(E) = S_{n,k} \tag 2$$
If $k\to \infty$ and $n$ is fixed, then the probability that an urn gets more than two balls turns negligible, and we can just count the events that have $n$ urns with $2$ balls and the rest empty. Summing the probabilities of these events we get
$$S_{n,k} \approx \binom{k}{n}\frac{(2n)!}{2^n} \tag 3$$
This approximation is similar to that of Mike Earnest, the graph displays $\log S_{n,k}$ for $n=18$ ("Ap1" is approx $(3)$, "Ap2" is Mike Earnest's)
It's also a lower bound, of course, because it omits some valid configurations.

Plugging in $(3)$ the Stirling approximation (first order) we get
$$S_{n,k} \approx (2 k n/e)^n$$
in agreement with your guess.
For both $n,k$ growing together, this approach (for Stirling numbers of the second kind) could be adapted.
Added: a really tight asymptotic, especially for large $k$ and $n$.
Let $m=2n$ and $\mu=m/k$, let $\lambda$ be the solution of $\lambda \tanh(\lambda) = \mu$ [*], and $\sigma^2=\lambda^2 + \mu(1-\mu)$
Then
$$S_{n,k} \approx \frac{ m! \, \cosh^k(\lambda)}{\lambda^{m}}
\sqrt{\frac{2}{ \pi k \sigma^2}} \tag 4$$
This approximation (lets call it "Ap3") is based on the same approach given in the link above (kind of Poissonization), I can provide details if anyone is interested.
It gives values that are practically indistinguishable from the exact values in the graph above.
Here's a table of values of $\log S_{n,k}$, for $n=18$.

and for $n=5$

[*] To numerically find $\lambda$, one can start with $\lambda_0=\max(\mu,\sqrt{\mu})$ and iterate $\lambda_{i+1}=\sqrt{\mu \lambda_i/\tanh(\lambda_i)}$. It converges quite fast.