Using argument similar to this answer, we can prove that the probability evaluates to $\frac{a^2-1}{3ab}$ for $2b\ge a>1$.
Let
$$x=2\sin\beta, y=2\sin(\alpha+\beta), z=2\sin\alpha$$
With
$$\alpha,\beta>0;\alpha+\beta<\pi$$
And consider
$$s=\frac{x}{z}=\frac{\sin\beta}{\sin\alpha},t=\frac{y}{z}=\frac{\sin(\alpha+\beta)}{\sin\alpha}$$
The mapping $(\alpha,\beta)\mapsto(s,t)$ is bijective onto the set of pairs $(s,t)\in\mathbb{R}^2$ satisfying
$$s,t>0;s+t>1>|s-t|$$
And the Jacobian determinant of the mapping is $st$, the condition $x,y,z$ be the sides of a triangle in random order and the vertices of the triangle be uniformly random on a circle has Lebesgue measure $\pi^2/2$.
Using the change of variable $(\alpha,\beta)\mapsto(s,t)$, we get
$$P\left(\frac{x^a+y^a}{(x+y)^a}<\frac{z^b}{(x+y)^b}\right)=\frac{2}{\pi^2}\int_{\Omega(a,b)}\frac{\mathrm ds\,\mathrm dt}{st}$$
With
$$\Omega(a,b)=\left\{(s,t)\in\mathbb{R}^2:s,t>0;s^a+t^a<(s+t)^{a-b}\right\}$$
Since beyond $s+t>1$, the region $s^a+t^a<(s+t)^{a-b}$ lies inside $|s-t|<1$, therefore we do the substitution:
$$m=s+t,n=s-t$$
And the probability becomes
$$\frac{2}{\pi^2}\cdot\color{red}2\cdot2\int_{\omega(a,b)}\frac{\mathrm dm\,\mathrm dn}{m^2-n^2}=\frac{8}{\pi^2}\int_{\omega(a,b)}\frac{\mathrm dm\,\mathrm dn}{m^2-n^2}$$
With
$$\omega(a,b)=\left\{(m,n)\in\mathbb{R}^2:m>1,\color{red}{0<n}<1,\left(\frac{m+n}{2}\right)^a+\left(\frac{m-n}{2}\right)^a<m^{a-b}\right\}$$
Let $m=u(n)$ be the unique positive real root of $\left(\frac{m+n}{2}\right)^a+\left(\frac{m-n}{2}\right)^a=m^{a-b}$, then $\left(\frac{m+n}{2}\right)^a+\left(\frac{m-n}{2}\right)^a<m^{a-b}\implies m<u(n)$, therefore by Fubini's theorem:
$$\frac{8}{\pi^2}\int_{\omega(a,b)}\frac{\mathrm dm\,\mathrm dn}{m^2-n^2}=\frac{8}{\pi^2}\int_{0}^{1}\int_{1}^{u(n)}\frac{\mathrm dm}{m^2-n^2}\,\mathrm dn=\frac{4}{\pi^2}\int_{0}^{1}\frac{\ln\frac{u(n)-n}{u(n)+n}-\ln\frac{1-n}{1+n}}{n}\,\mathrm dn$$
Let
$$n=uv,v\ge 0$$
$$\implies \left(\frac{u+uv}{2}\right)^a+\left(\frac{u-uv}{2}\right)^a=u^{a-b}$$
$$\implies u^{-b}=\left(\frac{1+v}{2}\right)^a+\left(\frac{1-v}{2}\right)^a$$
$$\implies u=\sqrt[-b]{\left(\frac{1+v}{2}\right)^a+\left(\frac{1-v}{2}\right)^a}$$
$$\implies n=v\sqrt[-b]{\left(\frac{1+v}{2}\right)^a+\left(\frac{1-v}{2}\right)^a}$$
When $v$ runs from $0$ to $1$, $n$ also runs from $0$ to $1$, and $u(n)$ stays above $1$, so we may do the substutition
$$n=v\sqrt[-b]{\left(\frac{1+v}{2}\right)^a+\left(\frac{1-v}{2}\right)^a}$$
$$\implies \frac{\mathrm dn}{n}=\frac{1}{v}-\frac{a}{b}\cdot\frac{(1+v)^{a-1}-(1-v)^{a-1}}{(1+v)^a+(1-v)^a}\,\mathrm dv$$
$$\implies \int_{0}^{1}\frac{\ln\frac{u(n)-n}{u(n)+n}-\ln\frac{1-n}{1+n}}{n}\,\mathrm dn=-\frac{a}{b}\int_{0}^{1}\frac{(1+v)^{a-1}-(1-v)^{a-1}}{(1+v)^a+(1-v)^a}\ln\frac{1-v}{1+v}\,\mathrm dv$$
Next, do the substitution
$$v=\tanh w\implies \mathrm dv=\frac{\mathrm dw}{\cosh^2 w}, \ln\frac{1-v}{1+v}=-2w,\\\frac{(1+v)^{a-1}-(1-v)^{a-1}}{(1+v)^a+(1-v)^a}=\frac{\cosh w\sinh((a-1)w)}{\cosh(aw)}$$
And the probability becomes
$$\frac{8a}{b\pi^2}\int_{0}^{\infty}\frac{w\sinh((a-1)w)}{\cosh w\cosh (aw)}\,\mathrm dw=\frac{8a}{b\pi^2}\int_{0}^{\infty}w(\tanh(aw)-\tanh w)\,\mathrm dw$$
Let
$$I(a)=\int_{0}^{\infty}w(\tanh(aw)-\tanh w)\,\mathrm dw$$
$$\implies I'(a)=\int_{0}^{\infty}\frac{w^2}{\cosh^2(aw)}\,\mathrm dw$$
Finally, do the substitution $z=aw$
$$\implies I'(a)=\int_{0}^{\infty}\frac{z^2}{a^3\cosh^2 z}\,\mathrm dz$$
It is proven in this answer, by rewriting the integrand then expand the denominator $\frac{1}{\left(1+e^{-2z}\right)^2}$, that
$$\int_{0}^{\infty}\frac{z^2}{\cosh^2 z}\,\mathrm dz=\frac{\pi^2}{12}$$
Therefore
$$I(a)=\int_{1}^{a}I'(t)\,\mathrm dt=\int_{1}^{a}\frac{\pi^2}{12t^3}\,\mathrm dt=\frac{\pi^2}{24}\left(1-\frac{1}{a^2}\right)$$
Connecting all the pieces together, we obtain
$$\boxed{P\left(\frac{x^a+y^a}{(x+y)^a}<\frac{z^b}{(x+y)^b}\right)=\frac{8a}{b\pi^2}\cdot\frac{\pi^2}{24}\left(1-\frac{1}{a^2}\right)=\frac{a^2-1}{3ab}}$$
