I notice two mistakes in your derivation thus far. It seems that your formula for $r(\phi)$ does not describe an ellipse. The correct expression in polar coordinates is
$$
r = \frac{ab}{\sqrt{a^2\sin^2(\phi)+b^2 \cos^2(\phi)}}
$$
Secondly, once you are in polar coordinates, the arc-length formula involves the square of $\frac{dr}{d\phi}$, not the second derivative $\frac{d^2 r}{d\phi^2}$:
$$
dl = \sqrt{r^2+\left(\frac{dr}{d\phi}\right)^2}
$$
Edit
It looks like you can put it in terms of a single ellptic integral, I think. Check out here: http://integraltec.com/math/math.php?f=ellipse.html, they have the full derivation from Cartesian, and the integral expression they get is the incomplete elliptic integral of the second kind. The difference is that their angle, $\theta$, is not a physical angle on the ellipse, but defined as $\sin(\theta)=x/a$. So you should be able to express your angle $\phi$ in terms of their angle $\theta$ by algebraic/trigonometric relationships and identities, then map $\phi$ to $\theta$ and substitute into the elliptic integral. This isn't a full solution, but gives you the flavor of how you can get an expression in terms of the ellipse parameters, elliptic integrals, and the radial angle to give you the arc length.
This is just an outline of an answer.
Call the integral to be calculated $I(a)$. First write
$$I(a) = 2\int_0^{\pi/2} \sqrt{\cos^2 k + a^2 \sin^2 k} \, dk$$
by symmetry. The difficulty here is that you can't just apply Taylor's formula for $a \to 0$, because values of $k$ near $\pi/2$ make a significant contribution to the integral and $\cos k$ is small there.
Make the substitution $u = \cot k$. Then
$$I(a) = 2 \int_0^{+\infty} \frac{\sqrt{u^2 + a^2}}{(u^2 + 1)^{3/2}} \, du.$$
Now divide the integral into three parts on the intervals $[0,a]$, $[a,1]$ and $[1,+\infty)$. Then
make the substitution $v = u/a$ on the first interval, $s = u^2$ on the second, and $w = 1/u$ on the third. We get
$$I(a) = 2a^2 \int_0^1 \sqrt{1+v^2} (1 + a^2 v^2)^{-3/2} \, dv + 2\int_0^1 (1 + w^2)^{-3/2} \sqrt{1 + a^2 w^2} \, dw \\+ \int_{a^2}^1 \frac{1}{(1+s)^{3/2}} \sqrt{1 + \frac{a^2}{s}} \, ds$$
Now the idea is to expand each integrand into a series using the Taylor series for $(1 + x)^{1/2}$ and $(1+x)^{-3/2}$, and then integrate term by term. This will be legitimate because of uniform convergence. In the first integral, expand the second factor as a function of $av$. Do the same in the second integral with respect to $aw$. The third is more complicated because $a^2$ appears as a bound, but you can expand the integrand into a double series with respect to $s$ and $a^2/s$. The resulting series is quite complicated.
However, if we only want an estimate at the level of $O(a^2)$, we can note that the square root in the last integral is $1 + a^2/2s$ to within $O(a^4/s^2)$, so the error in the integral will be at most $O(a^2)$. If we make this approximation, we find
$$I(a) = 2 - a^2 \ln a + O(a^2)$$
To evaluate the $a^2$ term is more difficult. The contribution from the first integral is $\sqrt{2} + \operatorname{arsinh}(1)$. The contribution from the second is $\operatorname{arsinh}(1) - \frac{1}{2}\sqrt{2}$. The $a^2$ term in the third integral consists of a $-a^2$ from the first term of the series, a term $a^2[-\operatorname{arsinh}(1)+ \frac{1}{\sqrt{2}} + \ln 2 - 1]$ in the second term of the series, and the remaining terms with coefficient $\sum_{n \geq 2} \frac{1}{n-1}\binom{1/2}{n}$. This last series is $f(1)$, where $f(x) = \sum_{n \geq 2} \frac{1}{n-1}\binom{1/2}{n}x^{n-1}$. We have $f(0) = 0$ and $f'(x) = \sum_{n \geq 2} \binom{1/2}{n}x^{n-2} = \frac{1}{x^2} (\sqrt{1 + x} - 1 - x/2)$, so
$$f(1) = \int_0^1 \frac{1}{x^2} (\sqrt{1 + x} - 1 - x/2) \, dx = \frac{3}{2} - \sqrt{2} + \ln 2 - \operatorname{arsinh}(1).$$
Taking everything into account, we get
$$I(a) = 2 - a^2\ln a + a^2 (2\ln 2 - 1/2) + O(a^4 \ln a).$$
Given how simple the result is, I'll bet there's a simpler way to find it.
EDIT: For $a=0.0001$, the true value of the integral is $2.000,000,100,966,347,688$. The approximation obtained by the formula is $2.000,000,100,966,347,331$.
Best Answer
The formula should be $$\color{blue}{s(x)=a E\left(\sin ^{-1}\left(\frac{x}{a}\right)|1-\frac{b^2}{a^2}\right)}$$ So, for your case $(a=20,b=5,x=\frac{73625}{10000}=\frac{589}{80})$, the result is $$20 E\left(\sin ^{-1}\left(\frac{589}{1600}\right)|\frac{15}{16}\right) \approx 7.37381$$ while $$20 E\left(\frac{589}{1600}|\frac{15}{16}\right) \approx 7.20786$$ So, a typo !