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).
Let me formalize the problem and explain 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) 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
- 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 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$
