Let us try not to use the (brilliant, yet a bit cumbersome) solution given 4 years ago, mentioned in the question comment by Random Variable and find a special function solution instead.
Note that what a "closed-form" solution refers to is subject to historical variations: elliptic functions did not belong to the set of closed-form primitives before they were developed in the mid-19th century. Similarly, Fox-Wright and Fox-H functions may not "look like" closed-form primitives for an undergraduate student and will look so at a higher level. I am choosing the latter convention. Note that implementation in standard CAS is also to be considered. Mathematica, for one, has the function $\texttt{FoxH}$.
Considering now
$$I_n = \int_0^{\frac{\pi}{2}} \frac{dx}{1+\sin^n(x)} = \sum_{m=0}^\infty (G_{m,n} - H_{m,n})$$
in which $G_{m,n} = \int_0^{\frac{\pi}{2}} \sin^{2mn}(x)\,dx, H_{m,n} = \int_0^{\frac{\pi}{2}} \sin^{(2m+1)n}(x)\,dx$
(Edit)
Note that the sums cannot be separated out (as in the alternate harmonic series), as they individually diverge:
$$ H_{m,n} \sim \sqrt{\frac{e\pi}{2(2m+1)n}}, \quad G_{m,n} \sim \sqrt{\frac{e}{\pi mn}}$$
By Gradshtein et al. [2010, eq. 3.621.3 p. 397],
$$G_{m,n} = \frac{(2mn-1)!!}{(2mn)!!} \frac{\pi}{2} = \frac{\Gamma(2mn)}{2^{2mn}\Gamma(mn)\Gamma(mn+1)}\pi,$$
for $m > 0$, with the convention $(-1)!!=1$. This formula may be simplified by way of Legendre's duplication formula
$$ \Gamma(mn)\, \Gamma(mn+\frac{1}{2}) = \Gamma(2mn) \; 2^{-2mn+1} \sqrt{\pi}$$
whence:
$$G_{m,n} = \frac{\sqrt{\pi}\Gamma\left(mn+\frac{1}{2}\right)}{2 \Gamma(mn+1)},$$
which has the advantage of being also true for $m=0$, and by eq. 3.621.4,
$$H_{m,k} = \frac{((2m+1)n-1)!!}{((2m+1)n)!!} \delta$$
\begin{equation}
H_{m,k} = \begin{cases}\frac{\sqrt{\pi} \Gamma((m+\frac{1}{2})n + \frac{1}{2})}{2 \Gamma((m+\frac{1}{2})n + 1)}\qquad \text{if $n$ is even}\\
2^{n-1}\frac{2^{2mn}\Gamma(mn+\frac{n-1}{2})^2}{\Gamma(2mn+n+1)} \qquad \qquad \qquad \text{if $n$ is odd}
\end{cases}
\end{equation}
for $m \geq 0$, in which $\delta= \frac{\pi}{2}$ if $n$ is even and $\delta=1$ if $n$ is odd.
As the sums of $(H_{m,n})_{m\geq 0}$ and $(G_{m,n})_{m\geq 0}$ individually diverge but converge if $(I_{m,n}=G_{m,n}-H_{m,n})_{m\geq 0}$ is considered, and the form of $H_{m,n}$ is very different from that of $G_{m,n}$ when $n$ is odd, we will be able to find a solution only for even values of $n$ using the special function method outlined below. Further developments should derive $I_n$ for odd values of $n$ from the preceding even values.
Considering therefore sums with even $n$, we can rewrite it as a Fox-Wright function (and not as a q-hypergeometric series as I first thought):
$$\sum_{m=0}^\infty I_{m,n}=\pi \, \left(\frac{1}{2} + \sum_{p=1}^\infty \frac{\Gamma(pn)}{2^{pn} \Gamma(\frac{p}{2}n)\Gamma(\frac{p}{2}n+1)}(-1)^p\right) = \,\pi\, {_2}\Psi_{2}\left[\begin{matrix}(0,n)&(1,1)\\ (0,\frac{n}{2})&(1,\frac{n}{2})\end{matrix}; -2^{-n}\right]$$
There is some basic algebra in the above formulae, which I still have not thoroughly checked (I will come back for edits if need be).
Hence the closed form expressed for $n$ even.
Note: Convergence conditions of the Fox-Wright function have not been checked, but my take is that the convergence of Fox-Wright functions is ensured by the existence of the integrals. For a closer view on convergence you may find it useful to check Mathai et al.'s The H function, Springer 2010, page 30.
To implement it using the Mathematica $\texttt{FoxH}$ function, simply use the fact that
$${_2}\Psi_{2}\left[\begin{matrix}(0,n)&(1,1)\\ (0,\frac{n}{2})&(1,\frac{n}{2})\end{matrix}; -2^{-n}\right] = H^{1,2}_{2,3}\left[\begin{matrix}(1,n)&(0,1)\\ (0,1) & (1,\frac{n}{2})&(0,\frac{n}{2})\end{matrix}; 2^{-n}\right]$$
