Hemisphere only. Only need to change limits of integration and leading weight of integral to solve for qarter sphere case.
$d^2=(r_1\cos\theta_1\sin\phi_1-r_2\cos\theta_2\sin\phi_2)^2+(r_1\sin\theta_1\sin \phi_1-r_2\sin\theta_2\sin\phi_2)^2+(r_1\cos\phi_1-r_2\cos\phi_2)^2$
$d^2=r_1^2\cos^2\theta_1\sin^2\phi_1+r_2^2\cos^2\phi_2\sin^2\phi_2^2-2r_1r_2\sin\phi_2\sin \phi_1\cos\theta_1\cos\theta_2+r_1^2\sin^2\theta_1\sin^2\phi_1+r_2^2\sin^2\theta_2\sin^2\phi_2-2r_1r_2\sin\theta_1\sin\theta_2\sin\phi_1\sin\phi_2+r_1^2\cos^2\phi_1+r_2^2\cos^2\phi_2-2r_1r_2\cos\phi_1\cos\phi_2$
$d^2=r_1^2+r_2^2-2r_1r_2\sin\phi_1\sin\phi_2\cos(\theta_1-\theta_2)-2r_1r_2\cos\phi_1\cos\phi_2$
$\frac{3}{2\pi R^3}\int_0^R\int_0^{2\pi}\int_0^{\pi/2}\sqrt{r_1^2+r_2^2-2r_1r_2\sin\phi_1\sin\phi_2\cos(\theta_1-\theta_2)-2r_1r_2\cos\phi_1\cos\phi_2} \ r_2^2\sin \phi_2 d\phi_2 d\theta_2dr_2$
WLOG $\theta_1=0.$
$\frac{3}{2\pi R^3}\int_0^R\int_0^{2\pi}\int_0^{\pi/2}\sqrt{r_1^2+r_2^2-2r_1r_2\sin\phi_1\sin\phi_2\cos(\theta_2)-2r_1r_2\cos\phi_1\cos\phi_2} \ r_2^2\sin \phi_2 d\phi_2 d\theta_2dr_2$
$r_1=\alpha R$. $r_2=\beta R. dr_1=Rd\alpha. dr_2=Rd\beta $
$F(\alpha, \phi_1)=\frac{3R}{2\pi }\int_0^1\int_0^{2\pi}\int_0^{\pi/2}\sqrt{\alpha^2+\beta^2-2\alpha\beta\sin\phi_1\sin\phi_2\cos\theta_2-2\alpha \beta\cos\phi_1\cos\phi_2}\beta^2\sin \phi_2d\phi_2d\theta_2d\beta$
$<d>=\frac{4}{\pi }\int_0^1\int_0^{\pi/2} F(\alpha,\phi_1) \alpha d\alpha d\phi_1$
I don't think the integrals have a closed form. Probably need to use elliptical integrals or numerical methods to evaluate.
Average distanced is $<d>=6CR/\pi^2$ where $C$ results from the integration.
Is interesting this is proportional to $1/\zeta(2)$, the probability that any two positive integers are relatively prime.