You’re most likely aware of this, but I’ll repeat it for the general reader to avoid spreading common misconceptions: As in the linked post, the term “random integers” is formalized by considering uniformly random positive integers up to some $n\in\mathbb N$ and taking the limit $n\to\infty$.
These numbers are always coprime exactly if they are never all divisible by the same prime; so the desired probability for them to always be coprime is the product over all primes $p$ of the probability that they are never all divisible by $p$. Thus we need to calculate the latter for general $p$.
Reduce the $u_{i,j}$ modulo $p$ and consider the $r$ vectors $x_j$ with $(x_j)_i=u_{i,j}$ in $\mathbb F_p^s$. Then we want the probability that $x_1+\sum_{j=2}^rx_jn_j$ is never the zero vector, that is, that the subspace spanned by $r-1$ vectors independently uniformly drawn from $\mathbb F_p^s$ doesn't contain the vector $-x_1$ also uniformly randomly drawn from $\mathbb F_p^s$. This is
$$
1-\sum_{d=0}^s\mathsf P(\text{the subspace has dimension $d$})p^{d-s}\;.
$$
I’ll first treat the case $r=2$ for general $s$, which includes your zero result for $s=2$ and your specific question for $s=3$. For $r=2$, the probability that the subspace has dimension $0$ is $p^{-s}$, and the probability that it has dimension $1$ is $1-p^{-s}$. Thus the desired probability is
\begin{eqnarray}
\prod_p\left(1-p^{-s}\cdot p^{0-s}-(1-p^{-s})\cdot p^{1-s}\right)
&=&
\prod_p\left(1-\frac1{p^{s-1}}+\frac1{p^{2s-1}}-\frac1{p^{2s}}\right)
\\
&=&
\prod_p\left(1-\frac1{p^{s-1}}\right)\prod_p\left(1+\frac{\frac1{p^{2s-1}}-\frac1{p^{2s}}}{1-\frac1{p^{s-1}}}\right)
\\
&=&
\prod_p\left(1-\frac1{p^{s-1}}\right)\prod_p\left(1+\frac{p-1}{p^{2s}-p^{s+1}}\right)
\\
&=&
\zeta(s-1)^{-1}\prod_p\left(1+\frac{p-1}{p^{2s}-p^{s+1}}\right)\;,
\end{eqnarray}
which is well approximated by $\zeta(s-1)^{-1}$. For $s=2$, the probability is zero since the zeta function has a pole at $1$, corresponding to the product $\prod_p\left(1-\frac1p\right)$ diverging to $0$.
For $s=3$, it is
$$
\zeta(2)^{-1}\prod_p\left(1+\frac{p-1}{p^6-p^4}\right)=\frac6{\pi^2}\prod_p\left(1+\frac1{p^4(p+1)}\right)\;.
$$
The product (which isn’t one of the named Euler products listed in Wikipedia) evaluates to approximately $1.02432$, so the answer to your specific question is about $62.3\%$.
I don’t know whether the calculation can be done in closed form for general $r$, but I’ll do it for $r=3$ to illustrate the principle.
We start again with $x_2$ spanning a subspace of dimension $0$ with probablity $p^{-s}$ and dimension $1$ with probability $1-p^{-s}$. If we now add another vector $x_3$, if the dimension was $0$ it remains $0$ with probability $p^{-s}$ and is increased to $1$ with probability $1-p^{-s}$, whereas if it was $1$ it remains $1$ with probability $p^{1-s}$ and is increased to $2$ with probability $1-p^{1-s}$. Thus, the desired probablity is
\begin{eqnarray}
&&
\prod_p\left(1-p^{-s}\cdot p^{-s}\cdot p^{0-s}-\left(p^{-s}\cdot(1-p^{-s})+(1-p^{-s})\cdot p^{1-s}\right)\cdot p^{1-s}-(1-p^{-s})(1-p^{1-s})\cdot p^{2-s}\right)
\\
&=&
\prod_p\left(1-\frac1{p^{s-2}}+\frac1{p^{2s-3}}-\frac1{p^{2s-1}}-\frac1{p^{3s-3}}+\frac1{p^{3s-2}}+\frac1{p^{3s-1}}-\frac1{p^{3s}}\right)
\\
&=&
\prod_p\left(1-\frac1{p^{s-2}}+(p^2-1)\left(\frac1{p^{2s-1}}-\frac1{p^{3s-1}}+\frac1{p^{3s}}\right)\right)
\\
&\approx&
\zeta(s-2)^{-1}\;.
\end{eqnarray}
Generally, for $r\ge s$ we get at least a term $\frac1p$ and thus the product diverges to $0$, whereas for $r\lt s$ the product is dominated by the term $\frac1{p^{s-r+1}}$ and well approximated by $\zeta(s-r+1)^{-1}$.
