Definite solution for the mean distance from an external point to the surface of a sphere. 
Sphere, radius $E$, is centred at point $O$ $[0,0,0]$.
External point $Q$ is at $[D,0,0]$. 
I can slice the sphere by making multiple planar cuts parallel to the $YZ$ plane to produce circular zones (quasi-discs) of equal infinitessimal width $dx$ and position $x_i$.  The curved surface of each zone has the same surface area.  The radius $R_i$ of any zone $i$ lies in the $YZ$ plane and has magnitude: $R_i = \sqrt{E^2 - x_i^2}  .$
The distance $L_i$ from point $Q$ to any point $P_i$ in the zone $i$ is given by:-
$$L_i = \sqrt{ (D-x_i)^2 + R_i^2} 
 = \sqrt{ D^2 -2Dx_i + x_i^2 + E^2 - x_i^2} 
= \sqrt{ D^2 -2Dx_i + E^2 } $$
$$L_i  = D \sqrt{ 1 -\frac{2x_i}{D} + \frac{E^2}{D^2} }. $$
Now I wish to integrate this expression over the range $-E\cdots+E$ and then divide by $2E$ to obtain the average value $\bar{L}$ thus
$$\bar{L} = \frac{D}{2E} \int_{-E}^{+E} \sqrt{ 1 -\frac{2x_i}{D} + \frac{E^2}{D^2} } dx.$$
I can make a Taylor Series approximation of the square root term using the standard formula $(1+x)^{0.5} = 1 + x/2 - x^2/8 + x^3/16 -\cdots$  
Applying this to the square root term I obtain
$$
\sqrt{ 1 -\frac{2x_i}{D} + \frac{E^2}{D^2} } = 1-x_i/D + (1/2)E^2/D^2 -(1/2)x_i^2/D^2 + (1/2)x_iE^2/D^3  + \cdots$$ where the subsequent terms on the RHS diminish in magnitude.  By dropping terms with odd powers of $x_i$ (because
they will go to zero when integrating over the range $-E \le x_i \le +E$ ) and consolidating and integrating I come up with the following approximate result:-
$$\bar{L} = \frac{D}{2E} \int_{-E}^{+E} \sqrt{ 1 -\frac{2x_i}{D} + \frac{E^2}{D^2} } dx
\approx
D\left(1 + \frac{E^2}{3D^2} + A\right) $$
where A is a small term whose expression depends on the number of terms evaluated in the Taylor series approximation.  From numerical modelling it appears plausible that the term $A$ vanishes if the Taylor Series is extended to infinite terms and thus the definite solution would be given by:-
$$\bar{L} = 
D\left(1 + \frac{E^2}{3D^2} \right)
= D + \frac{E^2}{3D}
$$
MY QUESTION
Is there a way to derive a definite solution to this problem?
 A: To complete the spheric coordinate computation, recall working with the unit sphere center origin, radius $1$ [this can be rescaled to fit the more general sphere radius]. And let the point outside the sphere be $P=(0,0,k)$ where $k>1$ ensures it is outside the sphere. The distance from any point $(x,y,z)$ to point $P$ is then $$f(x,y,z)=\sqrt{x^2+y^2+(k-z)^2}.$$ One then wants to integrate this over the sphere with parametrizion 
$$x=\sin t \cos \theta,\ y=\sin t \sin \theta, \ z=\cos \theta.$$
One must not forget to multiply $f(x,y,z)$ by the volume element $ \sin t$ in spherical coordinates as set up above. $f(x,y,z)$ becomes $\sqrt{k^2+1-2k \cos t}$ and after multiplying by $\sin t$ it becomes a simple integral via aubstitution, let $u$ be what's under the radical, etc. After integrating from $0$ to $\pi$ it becomes, if I'm not wrong, 
$$\frac{(k+1)^3-(k-1)^3}{3k}.$$
This could be further manipulated into $(6k^2+2)/(3k)$ but I thought it looked nicer in terms of the cube difference.
Added correction: In the above I have only integrated over $t \in [0,\pi]$ but must also integrate the resulting constant over $\theta \in [0,p\pi]$ to finish the integral. Then that integral must be divided by $4 \pi$ i.e. the surface area of the sphere of radius $1.$ The net result is to multiply my answer above by $1/2$ which makes it agree with steveOw's answer above. [Once my answer is adjusted to use the more general case of sphere radius $E$ and distance to sphere $D$ the two answers match.]
A: Following suggestions from coffeemath & user2566092 I found the following standard rule:-
$$
\int (ax + b) ^{p/n} = \frac{n}{(n+p)a} (ax + b) ^{1+p/n} +C
 $$
for $ p = \pm1, \pm2, \cdots p \ne -n .$
(Source: The Universal Encyclopedia of Mathematics, page 590; a similar equation is here).
Applying that to the expression for $\bar {L}$
$$\bar{L} = \frac{D}{2E} \int_{-E}^{+E} \sqrt{ 1 -\frac{2x_i}{D} + \frac{E^2}{D^2} } dx
=
 \frac{D}{2E} \int_{-E}^{+E} \left[ \left(\frac{-2}{D}\right)x_i  +    \left(1 + \frac{E^2}{D^2}\right)  \right]^{1/2} dx$$
whence $a=-2/D,b = 1 +(E^2/D^2), p = 1, n = 2$ and so
$$\bar{L} = 
 \frac{D}{2E} 
\frac{2}{3(-2/D)}
\left[ 
\left(
\left(\frac{-2}{D}\right)x_i  +    \left(1 + \frac{E^2}{D^2}\right)  \right)^{3/2}+constant\right]_{-E}^{+E} $$
then cancelling and dropping the constant leads to
$$\bar{L} = 
\frac{-D^2}{6E}
\left[ \left(
\frac{-2Dx_i+D^2+E^2}{D^2}   \right)^{3/2}\right]_{-E}^{+E} 
 = 
\frac{-D^2}{6E} \left( \frac{1}{D^2}\right)
^{3/2}
\left[ \left(
-2Dx_i+D^2+E^2  \right)^{3/2}\right]_{-E}^{+E} $$
expanding and shuffling minus signs
$$\bar{L} = 
\frac{1}{6DE}
\left[ 
\left(2DE+D^2+E^2  \right)^{3/2}
-
\left(-2DE+D^2+E^2  \right)^{3/2}
\right]
$$
this simplifies to
$$\bar{L} = 
\frac{1}{6DE}
\left[ 
\left( (D+E)^2  \right)^{3/2}
-
\left((D-E)^2 \right)^{3/2}
\right]
 = 
\frac{1}{6DE}
\left[ 
 (D+E)^3  -(D-E)^3 
\right]
$$
and then
$$\bar{L} = 
\frac{1}{6DE}
\left[ 
 6D^2E +2E^3 
\right]
=
D + \frac{E^2}{3D}
$$
and finally
$$\bar{L} = 
D \left( 1+ \frac{E^2}{3D^2} \right).
$$
