$\newcommand{\R}{\mathbb{R}}
\renewcommand{\P}{\operatorname{\mathsf P}}
\newcommand{\Ga}{\Gamma}
\newcommand{\de}{\delta}$
In view of the spherical symmetry of the distribution of the $l$-dimensional subspace, we can fix it to be, say, the span of the first $l$ vectors of the standard basis of $\R^n$ and, accordingly, let $v=:(Y_1,\dots,Y_n)$ be a random vector uniformly distributed on the unit sphere $S_{n-1}$.
So, the probability in question equals
\begin{equation}
p_\de:=\P\Big(\sum_1^l Y_i^2\le\de^2\Big);
\end{equation}
we are assuming that $\de\in(0,1)$.
Next, we may write
\begin{equation}
Y_i=X_i\Big/\sqrt{\sum_1^n X_i^2},
\end{equation}
where the $X_i$'s are iid $N(0,1)$, whence
\begin{equation}
p_\de=\P\Big(\frac{\sum_1^l X_i^2}{\sum_{l+1}^n X_i^2}\le\frac{\de^2}{1-\de^2}\Big)
=F_{l,n-l}\Big(\frac{\de^2}{1-\de^2}\frac{n-l}l\Big)
\sim c_{n,l}\de^l
\end{equation}
as $\de\downarrow0$,
where $F_{l,n-l}$ is the cdf of the $F$ distribution with $l,n-l$ degrees of freedom and
\begin{equation}
c_{n,l}:=\frac2{l\,B\left(l/2,(n-l)/2\right)}
=\frac{2\Ga(n/2)}{l\,\Ga(l/2)\Ga((n-l)/2)}.
\end{equation}
(Here we used the easily verifiable fact that, if the pdf $f$ of a nonnegative r.v. $X$ is such that for some real $c,p>0$ we have $f(x)\sim cx^{p-1}$ as $x\downarrow0$, then for the cdf $F$ of $X$ we have $F(x)\sim cx^p/p$ as $x\downarrow0$.)
