RULES ... don't just use divergent series to get an answer unless there is some additional work to show it is meaningful. Indeed, don't even use rearrangement of conditionally convergent series and expect to get something meaningful (again) unless there is some additional work.
AN EXAMPLE ... Fourier series sometimes diverge, but converge when summed (C,1). There are theorems explaining conditions when the (C,1) sum of a Fourier series for a function actually converges to that function.
Let me denote your LHS by $f(n,s)$. For fixed even $n$ I shall show that $f(n,s)-1\sim\zeta(s+1)-1$ as $s\to\infty$, that is,
$$\lim_{s\to\infty}\frac{f(n,s)-1}{\zeta(s+1)-1}=1.$$
This result nicely expresses your numerical observations, which show that the parts after the decimal point seem to be asymptotically the same.
On one hand, we have $\zeta(s+1)-1=2^{-s-1}+3^{-s-1}+\dots$. The terms after the second can be estimated from above by the integral $\int_2^\infty x^{-s-1}dx=\frac{2^{-s}}{s}$, so we see that $\zeta(s+1)-1\sim 2^{-s-1}$.
On the other hand, among pairs $(k,i)$ with $1\leq k\leq n,1\leq i\leq k$, the expression $\frac{\gcd(k,i)}{\operatorname{lcm}(k,i)}$ is equal to $1$ for exactly $n$ pairs $(k,k)$, and is equal to $2^{-1}$ for exactly $n/2$ pairs $(2k,k)$. All other terms, of which there are certainly fewer than $n^2$, are at most $3^{-1}$. Therefore we find
$$f(n,s)=\frac{1}{n}\left(n\cdot 1+\frac{n}{2}\cdot 2^{-s}+O(n^23^{-s})\right)=1+2^{-s-1}+o(2^{-s})$$
proving $f(n,s)-1\sim 2^{-s-1}$. It follows that $f(n,s)-1\sim\zeta(s+1)-1$, as we wanted.
Let me emphasize that in the above calculation it was crucial that $n$ was even. If $n$ is odd, then we instead only get $\frac{n-1}{2}$ pairs $(2k,k)$ and the asymptotics get slightly skewed - we then get $f(n,s)-1\sim\frac{n-1}{n}(\zeta(s+1)-1)$. For large $n$ the difference is however, pretty negligible.
Best Answer
The answer can be obtained with the following interpretation of the Ramanujan summation:
The function $R(x)=\zeta(z,x)+C\;$ satisfy $R(x) − R(x + 1) = x^{-z}$, $x>0$, $z\in \mathbb C$, $z\ne-1$, where $$ \zeta(z,x)=\sum_{k=0}^\infty\frac1{(k+x)^z} $$ is the Hurwitz zeta function (or its analytic continuation for $\Re z\le 1$.) The value of $C\;$ can be found with the help of the shift formula for the derivative $\frac{\partial}{\partial x} \zeta(z,x)=-z \zeta(z+1,x)\;$: $$ \int_1^2 R(x)dx= \int_1^2(\zeta(z,x)+C)dx= \int_1^2 \left( -\frac1{z-1}\frac{\partial}{\partial x} \zeta(z-1,x)+C\right)dx= $$ $$ C-\frac1{z-1}(\zeta(z-1,2)-\zeta(z-1,1))=C+\frac1{z-1}=0.$$ Hence $C=-\frac1{z-1}$. Also $\zeta(z,1)=\zeta(z)$ and we have $$ \sum_{n=1}^\infty n^{-z}=\zeta(z)-\frac1{z-1}(\mathfrak{R}),\quad z\ne1. $$ So Ramanujan summation transforms the Riemann zeta function into the regularized zeta function. It explains why the value $\gamma$ should be expected for the summation of the harmonic series.