Reynolds-averaged Navier-Stokes equations allow to split the description of a turbulent fluid into an averaged (typically laminar) flow on some length and/or time-scale and separate equations for the turbulent fluctuations. The resulting equations look like this $$\rho\bar{u}_j \frac{\partial \bar{u}_i }{\partial x_j} = \rho \bar{f}_i + \frac{\partial}{\partial x_j} \left[ - \bar{p}\delta_{ij} + \mu \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right) - \rho \overline{u_i^\prime u_j^\prime} \right ]$$ Here the bars denotes the averaged values, and $\rho \overline{u_i^\prime u_j^\prime}$ is called the Reynolds stress tensor characterizing the influence of the turbulent fluctuations on the mean flow.
To actually evaluate the Reynolds stress one usually passes to something called the Boussinesque hypothesis and that is that you can actually model the stress as a viscous stress tensor with a "turbulent viscosity" $\mu_t$ and an isotropic stress coming from "turbulent kinetic energy" $k$. That is, in Cartesian coordinates $$\rho \overline{u_i^\prime u_j^\prime} = \mu_t \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right) -\frac{2}{3} k \delta_{ij}$$ Then there is a number of models for how to compute the quantities $\mu_t$ and $k$. Just for illustration, one of them is the k-epsilon model, where two transport equations for variables $k,\epsilon$ are solved $$\frac{\partial (\rho k)}{\partial t}+ \frac {\partial (\rho k \bar{u}_i)}{\partial x_i}= f(k, \epsilon, \bar{u}_j, \partial \bar{u}_j/\partial x_k,...)$$ $$\frac{\partial (\rho \epsilon)}{\partial t}+ \frac {\partial (\rho \epsilon \bar{u}_i)}{\partial x_i}= g(k, \epsilon, \bar{u}_j, \partial \bar{u}_j/\partial x_k,...)$$ and turbulent viscosity is then determined as $\mu_t = \mu_t(k,\epsilon)$. Many other models exist.
Of course, in astrophysics we are talking about plasma dynamics, which is modeled by (radiative) compressible magneto-hydrodynamics. However, this set of equations can be Reynolds-averaged in very much the same way as the pure-fluid equations. The equations of models such as the k-epsilon model would have to be generalized by introducing the production of turbulent kinetic energy due to effects such as the magneto-rotational instability but otherwise the models should work in a similar fashion. Possibly, one would also have to include a model for the turbulent magnetic-field fluctuations in the Maxwell stress $\sim \overline{B_i B_j}$.
So now for my question: These Reynolds averaged models seem to have applications only in engineering, but I have never seen them applied in an astrophysical context. Why is this so?
I have instead seen a single, very special model, and that is the Shakura-Sunyaev prescription for turbulent viscosity in steady, thin accretion disks: $\mu_t = \alpha \rho \bar{p}$, where $\alpha$ is a constant. However, I do not see any other context than steady, thin disks where this kind of prescription can be useful. Does one perhaps use more sophisticated prescriptions in other astrophysical contexts such as the theory of stellar structure, intergalactic medium, or the solar wind?
