Let me formalize the problem and give some results. iv A number $n$ is $y$-smooth (or $y$-friable) if all its prime factors $p$ are less than or equal to $y$. Let us set $$\Psi(x,y) :=\# \{ n \le x: n \text{ is }y\text{-smooth}\},$$ which appeared in Sungjin Kim's comment. One can consider the quantity $$S_k(x,y) := \# \{ n_1,\ldots,n_k \le x: \gcd(n_1,n_2,\ldots,n_k)=1, \, n_i \text{ is }y\text{-smooth}\}.$$ The ratio $S_k(x,y)/\Psi(x,y)^k$ is exactly the probability that $n$ positive integers chosen uniformly at random from the set of $y$-smooth numbers up to $x$ have no non-trivial common factor. Equivalently, it is the probability that a vector of $y$-smooth numbers and $L^1$-norm $\le x$ is visible from the origin.
If $n \le x$ then $n$ is $x$-smooth, so $S_k(x,x)/\Psi(x,x)^k = S_k(x,x)/\lfloor x \rfloor^k$ is the probability that $n$ positive integers $\le x$ chosen uniformly at random have no common factor, and this converges to $1/\zeta(k)$. You seem to ask about $S_k(x,y)/\Psi(x,y)^k$ when $y$ is fixed.
We can use inclusion-exclusion to express $S_k(x,y)$ as $$\begin{align} S_k(x,y) &= \sum_{d \mid \prod_{p \le y}}\mu(d) \# \{ m_1,\ldots,m_k \le x/d: m_i \text{ is }y\text{-smooth}\}\\ &=\sum_{d\text{ is }y\text{-smooth}} \mu(d)\Psi(x/d,y)^k\end{align}$$ where $d$ stands for the common divisor of $k$ $y$-smooth numbers $n_1,\ldots,n_k$. So we need to understand $$ \sum_{d\text{ is }y\text{-smooth}} \mu(d)\frac{\Psi(x/d,y)^k}{\Psi(x,y)^k}.$$
A special case of a result of Ennola (1969, see discussion here) says that for fixed $y$, $\Psi(x,y) \sim C_y (\log x)^{\pi(y)}(1+O_y(1/\log x))$ where $\pi(y)$ is the number of primes up to $y$ and $C_y$ is an explicit positive constant.
Since $\log x \sim \log (x/d)$ for any fixed $d$, and since $\sum_{d\text{ is }y\text{-smooth}}\mu(d)=0$ for $y \ge 2$, this immediately implies that the above fraction is $O_y(1/\log x)$, so the probability tends to $0$. So most (in a limit sense) $k$-tuples of $y$-smooth integers have a common factor, as observed numerically by Stefan Kohl and Aaron Meyerowitz.
Below I discuss what happens if $y$ is allowed to grow with $x$. First of all, in this case it is often useful to introduce a truncation parameter $T>0$ into the above inclusion-exclusion formula, and obtain via the union bound the estimate $$ (\star)\, S_k(x, y) = \sum_{d \text{ is }T\text{-smooth}} \mu(d) \Psi(x/d,y)^k + O\left( \sum_{y\ge p > T} \Psi(x/p,y)^k \right).$$ The ratio $\Psi(x/d,y)/\Psi(x,y)$ (`local density') was studied extensively in the literature, namely in
- Adolf Hildebrand and Gérald Tenenbaum, "On integers free of large prime factors". Trans. Amer. Math. Soc. 296 (1986), no. 1, 265–290.
- Aleksandar Ivić and Gérald Tenenbaum, "Local densities over integers free of large prime factors". Quart. J. Math. Oxford Ser. (2) 37 (1986), no. 148, 401–417.
- Régis de la Bretèche and Gérald Tenenbaum, "Propriétés statistiques des entiers friables", Ramanujan J. 9 (2005), no. 1-2, 139–202.
- See also Theorem III.5.22 in Tenenbaum's book (3rd edition of English version).
These papers prove either asymptotic results or bounds on $\Psi(x/d,y)/\Psi(x,y)$ in terms of $d^{-\alpha}$ where $\alpha\in (0,1)$ is a certain 'saddle point' defined in terms of $x$ and $y$. Below is a quick lemma one can prove using the paper of Ivić and Tenenbaum, which is quite powerful when $y$ is not too small.
Lemma: Fix $k \ge 2$. Suppose $x \ge y \ge (\log x)^{1 + \frac{1}{k-1}+\varepsilon}$ for some $\varepsilon>0$. Then $$(\star \star)\, \frac{S_k(x,y)}{\Psi(x,y)^k} \sim \prod_{p}\left(1-p^{-k\beta(x,y)}\right)$$ as $x \to \infty$, where $\beta:=1-\frac{\xi(u)}{\log y}$, $u:=\log x / \log y$ and $\xi(u)\sim \log u$ is defined via $e^{\xi(u)}-1=\xi(u)$. In particular, the infinite product in $(\star \star)$ is bounded away from $0$.
Remark: When $y=x$ we have $\beta = 1$ and the limit is $1/\zeta(k)$ exactly.
Proof: Lemma 2 of the Ivić--Tenenbaum paper says that, for the above $\beta$, $$\Psi(x/d,y) = \Psi(x,y)d^{-\beta}\left(1 + O_{\varepsilon}\left( \frac{\log d}{\log x} + \frac{\log \log y}{\log y}\right)\right).$$ when $1 \le d \le y$ and $y \ge (\log x)^{1+\varepsilon}$. Lemma 3 says that $$\Psi(x/d,y) \ll \Psi(x,y) d^{-\beta + \frac{c}{\log y}}$$ uniformly for $1 \le d \le x$ and $y \ge (\log x)^{1+\varepsilon}$, and $c$ is an absolute positive constant. We use $(\star)$ with $T=(\log y)/2$ and these two estimates to obtain $$\frac{S_k(x, y)}{\Psi(x,y)^k} = (1+o(1))\prod_{p \le \log y/2} \left(1-p^{-k\beta}\right) + O\left(\sum_{y \ge p > \log y/2} p^{-\beta k}\right).$$ The assumption $y \ge (\log x)^{1+1/(k-1)+\varepsilon}$ implies $\beta k \ge 1+c\varepsilon$, and the lemma follows. $\blacksquare$