I have not proven the following to work, but I've done similar calculations when debugging software in my 2009-10 work with Matt Hastings. If these seem fast enough, it should not be too much work to validate them.

(1) I would take small, random perturbations of the identity matrix and then do four or five iterations of Newtons method to get the polar part of the resulting matrix. I offer this suggestion under the assumption that you can tolerate some samples falling outside the delta ball.

By Newton's method I mean replacing $X$ by $\frac{1}{2}\left( X + (X^{-1})^\mathrm{T}\right) $ where $T$ is the transpose. This is a simplification of what is suggested by N. J. Higham for computing the polar part of an invertible matrix.

The number of iterations needed depends on delta. I have no idea how big your matrices are. If that are bigger than 200 by 200 you might need to use a package that has a well parallelized matrix inversion routine. If you are looking at smaller matrices, you might just go ahead and apply 10 iterations.

(2) If you cannot abide by matrices outside the delta ball, then generate a random diagonal orthogonal $D$ that is close to $I$ and a random orthogonal $W$ and multiply to $U = WDW^\mathrm{T}$ as the desired random orthogonal, now known to be in the delta ball. I am assuming the operator norm, which is costly to compute by the default in pure math, I think.

You need to generate the random unitary by taking polar decomposition of a random matrix. Now you need to read up on Newton's method for for computing the polar part, as you may deal with badly conditioned matrices.

I am not sure what distribution you want for the diagonal orthogonal matrix. So I am hoping (1) will work for you.