Any way I can evaluate this integral?

This is not an answer, but a hint how to possibly proceed. It is a pity to pass on some nice Wolfram Language abilities. Let me know if I should remove it.

Idea is: I think YOU should be the first to derive this integral - if it is possible. Consider only rational numbers. If this integral behaves smoothly - which it obviously does:

DiscretePlot[NIntegrate[1/(Sin[s/2]^b + 1), {s, 0, Pi}], {b, 0, 3, .1}] 

we can consider integration only with rational exponents because they can fill between irrational and transcendental exponents pretty well to represent smooth behaviour. Then for some simple $\beta$ you'll get all sorts of formulas - no consistency:

Integrate[1/(Sin[s/2]^(1/2) + 1), {s, 0, Pi}] 

$2 \sqrt{2} E\left(\frac{1}{2}\right)-2$

Integrate[1/(Sin[s/3]^(1/3) + 1), {s, 0, Pi}] 

$\frac{\sqrt{\pi } \left(\left(88 (-1)^{3/8} \sqrt[6]{4-2 \sqrt{2}} \left(10 \left(58+41 \sqrt{2}\right) \, _2F_1\left(\frac{1}{6},\frac{1}{3};\frac{7}{6};17-12 \sqrt{2}\right)-\left(24+17 \sqrt{2}\right) \, _2F_1\left(\frac{1}{3},\frac{5}{12};\frac{17}{12};17-12 \sqrt{2}\right)\right)+(20+20 i) \sqrt[6]{2} \left(22 \left((17+7 i)+(12+5 i) \sqrt{2}\right)-7 \sqrt[6]{2} \sqrt[3]{3 \sqrt{2}-4} \left((7+3 i)+(5+2 i) \sqrt{2}\right) \, _2F_1\left(\frac{1}{3},\frac{11}{12};\frac{23}{12};17-12 \sqrt{2}\right)\right)\right) \Gamma \left(-\frac{2}{3}\right)+55 \left((41+99 i)+(29+70 i) \sqrt{2}\right) \left(\Gamma \left(\frac{1}{6}\right)^2+2 \Gamma \left(-\frac{1}{3}\right) \Gamma \left(\frac{2}{3}\right)\right)\right)}{55 \left((41+99 i)+(29+70 i) \sqrt{2}\right) \Gamma \left(-\frac{1}{3}\right) \Gamma \left(\frac{1}{6}\right)}$

which is pretty terrifying. But I think these are all simplifications form the Meijer G function as can be seen from more general $\beta$:

Integrate[1/(Sin[s/2]^(5/7) + 1), {s, 0, Pi}] 

Integrate[1/(Sin[s/2]^(7/5) + 1), {s, 0, Pi}] 

Try it with other fractions. So if you can figure out how arguments of Meijer G function are filled in you can probably publish a paper ;-)

Here is a series-expansion way of doing your integral.

Expand the integrand as

Sum[(-1)^n Sin[θ/2]^(n β), {n, 0, ∞}] (* 1/(1+Sin[θ/2]^β) *) 

Integrate each term in the sum using

Integrate[Sin[θ/2]^(n β), {θ, 0, π}] (* ConditionalExpression[(Sqrt[π] Gamma[1/2 (1 + n β)])/Gamma[1 + (n β)/2], Re[n β] > -1] *) 

Fold it all together to obtain

Sum[(-1)^n (Sqrt[π] Gamma[1/2 (1 + n β)])/Gamma[1 + (n β)/2], {n, 0, ∞}] 

This gives the same numerical result as direct evaluation using NIntegrate.

This is a way which will give you a the final solution in terms of $\beta$. The idea is generate a Table of solution for different $\beta$ value and then find a fitting function. I choose here $-10<\beta<10$.

data = Table[{\[Beta],NIntegrate[1/(Sin[\[Theta]/2]^\[Beta] + 1), {\[Theta],0,\[Pi]}]}, {\[Beta], -10, 10, .5}]; ListPlot[data] sol1 = Interpolation[data] ;(*Interpolation*) Plot[sol1[x], {x, -10, 10}, Epilog -> Point[data],PlotLabel->"sol1"] sol2[x_] = Fit[data, Table[x^n, {n, 0, 20}], x]; (*Polynomial fit*) Plot[sol2[x],{x,-10,10},Epilog->Point[data],PlotLabel->"sol2"] 

So both $\rm sol1[\beta]$ and $\rm sol2[\beta]$ are your solution which you can see from the plots. Regarding $\rm sol2[\beta]$ you better use higher order polynomial (I use here 20) to get better result.

Following the suggestion by @george2079, you can also guess an approximate function for your solution with some unknown parameter and get an approximate (close enough) solution by fitting it. Here I show it for a+b ArcTan[c x].

guess[x_] := a + b ArcTan[c x] A = FindFit[data, guess[x], {a, b, c}, x]; sol3 = guess[x] /. A Plot[sol3, {x, -10, 10}, Epilog -> Point[data], PlotLabel -> "sol3"] 

and your final solution is $\rm 1.5708 + 0.776552 Tan^{-1}(0.554812 x)$. As suggested by @george2079, I think the 1.5708 will converge to $\pi$ if you consider more points. Anyway you may not be lucky enough every time to get such a nice guess.