  Turbulence models attempt to account for unresolved dynamics and diffusion in
hydrodynamical simulations. We develop a common framework for two-equation
Reynolds-Averaged Navier-Stokes (RANS) turbulence models, and we implement six
models in the Athena code. We verify each implementation with the standard
subsonic mixing layer, although the level of agreement depends on the
definition of the mixing layer width. We then test the validity of each model
into the supersonic regime, showing that compressibility corrections can
improve agreement with experiment. For models with buoyancy effects, we also
verify our implementation via the growth of the Rayleigh-Taylor instability in
a stratified medium. The models are then applied to the ubiquitous
astrophysical shock-cloud interaction in three dimensions. We focus on the
mixing of shock and cloud material, comparing results from turbulence models to
high-resolution simulations (up to 200 cells per cloud radius) and
ensemble-averaged simulations. We find that the turbulence models lead to
increased spreading and mixing of the cloud, although no two models predict the
same result. Increased mixing is also observed in inviscid simulations at
resolutions greater than 100 cells per radius, which suggests that the
turbulent mixing begins to be resolved.
