This question, at least for $k=2$, was studied by Gunderson, Coppersmith and Granville. Surprisingly, the motivation has to do with the first case of Fermat's Last Theorem, but I won't elaborate on that. See specifically Granville's paper "On positive integers $\le x$ with prime factors $\le t \log x$" (Number Theory and Applications (ed R.A. Mollin) (Kluwer, NATO ASI, 1989), pages 403-422, [PDF link][1]).
Let me formalize the problem and explain some some results.
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][2]) 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.
-----
In Granville's paper linked above, he studies the quantity $\Psi(x,x',y)$, which counts pairs $(a,b)$ of coprime $y$-smooth integers with $a \le x$ and $b\le x'$. For $x'=x$, $\Psi(x,x,y)$ is precisely $S_2(x,y)$. In pages 8--9 he quickly derives the following results, and the arguments should extend to $S_k(x,y)$:
- For $2 \le y \le (\log x)^{1/2}$, we have $\Psi(x,x,y) \sim \binom{2\pi(y)}{\pi(y)} \Psi(x,y)$ as $x \to \infty$. Note that this implies that for fixed $y \ge 2$, the probability $S_2(x,y)/\Psi(x,y)^2$ decays like $1/\Psi(x,y) \asymp (\log x)^{-\pi(y)}$.
- For $x \ge y \ge 2$ with $y/\log x \to 0$ as $ x \to \infty$, we have $\Psi(x,x,y) = \Psi(x,y)^{1+o(1)}$.
- For $x \ge y \ge 2$ with $y/\log x \to \infty$ we have $\Psi(x,x,y)=\Psi(x,y)^{2+o(1)}$ as $x \to \infty$.
- For $x \ge y \ge (\log x)^{2+\varepsilon}$ we have $\Psi(x,x,y) \sim \Psi(x,y)^2/\zeta(2\alpha)$ as $x \to \infty$ where $\alpha$ is the saddle point associated with $x$ and $y$. (I'll include a proof below.)
This leaves out the range $y \asymp \log x$, which is the main focus of Granville's paper. In Theorem 3 he proves an asymptotic formula in this remaining case.
----
Below I explain how one treats the case where $y$ grows sufficiently fast with $x$ (namely, $y \ge (\log x)^{1+1/(k-1)+\varepsilon}$, but I'll start with a general discussion.
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 the strategy hinted by Granville.
----
**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 \zeta(k)^{-1}$$
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 \zeta(k\beta(x,y))^{-1}$$
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://dms.umontreal.ca/~andrew/PDF/Smooth1.pdf
[2]: https://mathoverflow.net/questions/420881/reference-request-numbers-composed-of-given-primes