That’s the way to do it. If you had a closed form of your posterior to compute the covariance from then you wouldn’t need to sample!
There was another post yesterday where someone was modeling the covariance directly as a parameter (using LKJ priors and such). The covariance in this case was a hyperparameter of some random effects, and the covariance was something of particular interest for the application. It’s pretty well written I think: Reparameterizing to avoid low E-BFMI warning - #5 by mdonoghoe .
To go beyond that your best bet would be looking through BDA3: http://www.stat.columbia.edu/~gelman/book/