$\newcommand{\ep}{\varepsilon}
\newcommand{\Ga}{\Gamma}
\newcommand{\de}{\delta}$
For $\ep\in(0,1)$, let
\begin{equation*}
P_{m,n}:=P\Big(\bigcap_{1\le i<j\le m}\{|v_i\cdot v_j|\le\ep\}\Big)
=1-Q_{m,n},
\end{equation*}
where
\begin{equation*}
Q_{m,n}:=P\Big(\bigcup_{1\le i<j\le m}\{|v_i\cdot v_j|>\ep\}\Big).
\end{equation*}
By Bonferroni inequalities,
\begin{equation*}
Mp\ge Q_{m,n}\ge Mp-R/2,
\end{equation*}
where
\begin{equation*}
p:=P(|v_1\cdot v_2|>\ep),
\end{equation*}
\begin{equation*}
M:=m(m-1)/2,
\end{equation*}
\begin{equation*}
R:=\sum_{1\le i<j\le m}\;\sum_{1\le k<l\le m,\,(k,l)\ne(i,j)}
P(|v_i\cdot v_j|>\ep,|v_k\cdot v_l|>\ep).
\end{equation*}
If $\{i,j\}\cap\{k,l\}=\emptyset$, then $P(|v_i\cdot v_j|>\ep,|v_k\cdot v_l|>\ep)=P(|v_i\cdot v_j|>\ep)\,P(|v_k\cdot v_l|>\ep)=p^2$.
If $\{i,j\}\cap\{k,l\}\ne\emptyset$ but $\{i,j\}\ne\{k,l\}$, then, using the iid condition on the $u_i$'s and the spherical symmetry, for (say) the unit vector $e_1$ of the standard orthonormal basis of $\mathbb R^n$, we have
\begin{multline}
P(|v_i\cdot v_j|>\ep,|v_k\cdot v_l|>\ep)
=P(|v_1\cdot v_2|>\ep,|v_1\cdot v_3|>\ep) \\
=P(|e_1\cdot v_2|>\ep,|e_1\cdot v_3|>\ep)
=P(|e_1\cdot v_2|>\ep)\,P(|e_1\cdot v_3|>\ep)=p^2.
\end{multline}
So, $P(|v_i\cdot v_j|>\ep,|v_k\cdot v_l|>\ep)=p^2$ for any $i,j,k,l$ such that $1\le i<j\le m,\,1\le k<l\le m,\,(k,l)\ne(i,j)$. So,
\begin{equation*}
R=M(M-1)p^2.
\end{equation*}
Next,
\begin{equation*}
p=P(|e_1\cdot v_1|>\ep)=K_nI_n,
\end{equation*}
where, with $n\to\infty$,
\begin{equation*}
K_n:=\frac{\Ga(n/2)}{\Ga(1/2)\Ga((n-1)/2)(n-1)^{1/2}}\to1/\sqrt\pi,
\end{equation*}
\begin{equation*}
I_n:=(n-1)^{1/2}\int_{\sqrt c}^\infty(1+t^2)^{-n/2}\,dt=e^{-nc/(2+o(1))},
\end{equation*}
\begin{equation*}
c:=\frac{\ep^2}{1-\ep^2}.
\end{equation*}
Collecting the pieces, we see that
\begin{equation}
Q_{m,n}=Me^{-nc/(2+o(1))}.
\end{equation}
Setting now $\de=Q_{m,n}\to0$, we find the asymptotics of the needed $n$:
\begin{equation}
n\sim2\frac{1-\ep^2}{\ep^2}\,\ln\frac{m(m-1)}{2\de}.
\end{equation}
