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][1]) 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
1. Adolf Hildebrand and Gérald Tenenbaum, "On integers free of large prime factors". Trans. Amer. Math. Soc. 296 (1986), no. 1, 265–290.
2. 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.
3. 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.
4. 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.
----
**Corollary**: Fix $k \ge 2$. Suppose $x \ge y$ and that $y$ grows with $x$ faster than any power of $\log x$ (i.e. $\log y /\log \log x \to \infty$). Then
$$\frac{S_k(x,y)}{\Psi(x,y)^k} \sim 1/\zeta(k)$$
as $x \to \infty$.
This is immediate from
**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$.
**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$ tending to infinity sufficiently slow (we want $d \le y$ and $\log d/\log x \to 0$. Since $d$ is at most $\prod_{p \le T}p=e^{T(1+o(1))}$ we can take, say, $T=\log y/\log \log y$), and these two estimates to obtain
$$\begin{align}
\frac{S_k(x, y)}{\Psi(x,y)^k} &= \prod_{p \le T} \left(1-p^{-k\beta}\right) + O\left(\sum_{y \ge p > T} p^{-\beta k}\right)\\
&+O_{\varepsilon}\left( \left( \frac{T}{\log x} + \frac{\log \log y}{\log y}\right) \prod_{p \le T} \left(1+p^{-k\beta}\right)\right).
\end{align}$$
The assumption $y \ge (\log x)^{1+1/(k-1)+\varepsilon}$ implies $\beta k \ge 1+c\varepsilon$, and the lemma follows. $\blacksquare$
[1]: https://mathoverflow.net/questions/420881/reference-request-numbers-composed-of-given-primes