About the limit $\lim_{n \to +\infty} \frac{1}{n^2} \sum_{1 \le a,b \le n} \frac{1}{ \mathrm{gcd} (a,b)} $

There is no need to involve probabilities, just a bit of analytic number theory suffices. I don't know how familiar you are with $L$-functions, so I'll leave some details open for now. If any part needs clarification, let me know.

For every positive integer $n$ define $$f(n):=\sum_{a,b=1}^n\frac{1}{\gcd(a,b)} \qquad\text{ and }\qquad g(n):=\sum_{d\mid n}\frac{\varphi(\tfrac{n}{d})}{d}$$ so that for all $n>1$ we have \begin{eqnarray*} f(n)-f(n-1) &=&-\frac1n+2\sum_{c=1}^n\frac{1}{\gcd(c,n)}\\ &=&-\frac1n+2\sum_{d\mid n}\frac{\varphi(\tfrac{n}{d})}{d}\\ &=&-\frac1n+2g(n). \end{eqnarray*} Because $\lim_{n\to\infty}\frac{1}{n^2}\sum_{k=1}^n-\frac1k=0$, this already shows that $$\lim_{n\to\infty}\frac{1}{n^2}f(n)=2\lim_{n\to\infty}\frac{1}{n^2}\sum_{k=1}^ng(n).$$ Note that $g$ is a multiplicative function, and in fact $g=\iota\ast\varphi$ where $\iota(n):=\tfrac1n$ and $\ast$ denotes the Dirichlet convolution of multiplicative functions. This immediately shows that the $L$-function associated to $g$ converges for all $s\in\Bbb{C}$ with $\operatorname{Re}(s)>2$ and that for all such $s\in\Bbb{C}$ we have $$L_g(s)=L_{\iota}(s)L_{\varphi}(s)=\zeta(s+1)L_{\varphi}(s).$$ It is easy to see that $L_{\iota}(s)=\zeta(s+1)$ for these $s$, and it is a basic exercise to show that $$L_{\varphi}(s)=\frac{\zeta(s-1)}{\zeta(s)},$$ whenever $\operatorname{Re}(s)>2$. Then by a Tauberian theorem, for example Wiener-Ikehara, we have \begin{eqnarray*} \lim_{n\to\infty}\frac{1}{n^2}f(n) &=&2\operatorname{Res}(L_g,2)\\ &=&2\cdot\operatorname{Res}(\zeta(s+1),2)\cdot\operatorname{Res}(L_{\varphi},2)\\ &=&2\cdot\zeta(3)\cdot\frac{1}{2\zeta(2)}\\ &=&\frac{\zeta(3)}{\zeta(2)}. \end{eqnarray*}

Your doubts are justified: $\sum_{d=1}^n \ \frac{\zeta(2)^{-1}}{d^2}$ is just a sum of limiting probabilities, so it should (and does) yield $1$ in the limit. The sum you want to compute is of values of $\frac{1}{\gcd(a,b)}$, so you need to weight the possible values of that expression by their probabilities. This gives $$ \sum_{d=1}^n \ \frac{\zeta(2)^{-1}}{d^2}\cdot\frac{1}{d}\to\frac{\zeta(3)}{\zeta(2)}\approx 0.730763 $$ A numerical experimentat (with $n=1000$) gives a result close to $0.731$, so this seems likely to be correct.