Actually the above is the formula for the exact field within the toroidal winding if we neglect things like finite wire diameter and circumferential current. You can see this because reflection through a plane containing the principal axis of symmetry of the toroid doesn't change the current distribution, so it doesn't change the magnetic field either. However such a mirror reflection reverses the in-plane components of $\vec B$ and leaves the out of plane components unchanged because $\vec B$ is an axial vector field. Thus only the out of plane, i.e. azimuthal component of $\vec B$ can be nonzero. Then rotational symmetry tells us that $B_{\phi}$ doesn't depend on the azimuthal angle $\phi$ so Ampere's law gets us exactly the $\vec B$ field.
To get flux through a loop of the coil we might use a formula like $\Phi_1=\frac{\mu_0NIA}{2\pi}\left(\frac1r\right)_{\text{average}}$ where $A$ is the cross-sectional area of the coil. For a toroid with circular cross section of radius $b$ centered a distance $a$ from the principal axis of symmetry, we have
$$\Phi_1=\int_0^{2\pi}\int_0^b\frac{\mu_0NI}{2\pi(a-r\cos\theta)}r\,dr\,d\theta=\frac{\mu_0NIb^2}{a+\sqrt{a^2-b^2}}=\frac{\mu_0NI\pi b^2}{2\pi}\left(\frac1r\right)_{\text{average}}$$
So for this geometry we get
$$\left(\frac1r\right)_{\text{average}}=\frac2{a+\sqrt{a^2-b^2}}$$
For a square cross section with side $2b$ centered a distance $a$ from the principal axis of symmetry we get
$$\Phi_1=\int_{-b}^b\int_{-b}^b\frac{\mu_0NI}{2\pi(a-x)}dx\,dy=\frac{2b\mu_0NI}{2\pi}\ln\left(\frac{a+b}{a-b}\right)=\frac{4b^2\mu_0NI}{2\pi}\left(\frac1r\right)_{\text{average}}$$
So for this geometry
$$\left(\frac1r\right)_{\text{average}}=\frac1{2b}\ln\left(\frac{a+b}{a-b}\right)$$
In either case the answer approaches $1/a$ as $b\rightarrow0$ but the exact form depends on coil geometry.