I arXived a paper on this topic recently https://arxiv.org/pdf/1906.07684.pdf
Essentially, to simulate from the distribution of a random orthogonal matrix Q, we define a distribution for a real random matrix X such that the orthogonal component of the polar decomposition of X is equal in distribution to Q. This can be viewed as a form of parameter expansion. For a particular choice of distribution for the expanded parameters, this approach looks a lot like what Stan does to simulate from the unit sphere. See Equation 4.
*I should add that there is a typo directly following Equation 4. The argument of the density should be Q, not Q_X.