$\int_{0}^{1}\int_{0}^{1}\frac{1}{\sqrt{x^2+y^2}}\,dx\,dy=$?

Here is another way to perform the integral without using polar coordinates.

Notice the integrand and the domain to integrate are both symmetric with respect to $x \leftrightarrow y$.
This means we can split the integral into two equal pieces. One piece for $x \le y$ and another piece for $x \ge y$. On one of the piece, say $x \le y$, introduce variable $t$ such that $x = yt$, we have:

$$ \int_0^1\int_0^1 \frac{1}{\sqrt{x^2+y^2}} dx dy = 2 \int_0^1 \int_0^y \frac{1}{\sqrt{x^2+y^2}} dx dy = 2 \int_0^1 \left(\int_0^1 \frac{1}{y\sqrt{1+t^2}} ydt\right) dy\\ = 2\int_0^1\frac{1}{\sqrt{1+t^2}}dt = 2\sinh^{-1}(1) = 2\log(\sqrt{1^2+1}+1) = 2\log(\sqrt{2}+1)$$

Your substitution works. With that, you should get: $$\int_0^1 \int_0^1 \frac{1}{\sqrt{x^2+y^2}}\,dx\,dy=2\int_0^{\pi/4}\int_0^{1/\cos\theta} dr\,d\theta =2\int_0^{\pi/4} \frac{d\theta}{\cos\theta}$$ Since: $$\int_0^{\pi/4} \frac{d\theta}{\cos\theta}=\int_0^{\pi/4} \sec\theta\, d\theta=\left(\ln|\sec \theta+\tan\theta|\right|_0^{\pi/4}=\ln(\sqrt{2}+1)$$ Hence, the final result is: $$\boxed{2\ln(\sqrt{2}+1)}$$

---

To get the bounds for $r$, consider the figure below:

From $\Delta OAC$, $$\cos\theta=\frac{OA}{OC}=\frac{1}{OC} \Rightarrow OC=\frac{1}{\cos\theta}$$ Hence, $r$ varies from $0$ to $1/\cos\theta$.

you made a good guess, changing to polar is the right thing to do $x=r\cos\theta$ and $y=r\sin\theta$ does the trick. Also probably you know that the Jacobian is $r$ do your integral becomes $$\iint \frac{1}{r}r dr d\theta$$ Since the region is a square, the bounds for theta is $0\to \pi/2$. The bounds for $r$ is varying.

So, we should split up the integral as $$\iint = \int_0^{\pi/4}\int_0^{r_1} + \int_{\pi/4}^{\pi/2}\int_0^{r_2}$$

$r_1 = \sqrt{1^2 + y^2} = \sqrt{1^2 + y^2} = \sqrt{1 + r_1 ^2 \sin^2 \theta}$ solve for $r_1$ from here ( you will get $r_1 = 1/\cos\theta $ ) and do the same for $r_2$ on the region $\mathrm{II}$