Athena Applications Page 
The most uptodate (but not necessarily complete) source of papers that describe applications that have used Athena can be found by browsing the list of papers that cite the ApJS method paper.
Images from, and links to, some of these papers are given below.

Three Dimensional Compressible Hydrodynamic Simulations of Vortices in Disks by Y. Shen, J. Stone, & T. Gardiner, ApJ, 653, 513 (2006). The image shows the component of the specific vorticity in the vertical direction from a threedimensional hydrodynamic simulation of the evolution of an initially random vorticity field in a Keplerian shear flow. In 3D the vorticity and kinetic enerrgy quickly die away. 

Effect of the Coriolis Force on the Hydrodynamics of Colliding Wind Binaries by N. Lemaster, J. Stone, & T. Gardiner, ApJ 662, 582 (2007). Shock structure in the colliding winds from two identical stars on circular orbits, with the orbital velocity 40% of the stellar wind velocity. Small spheres mark the location of the stars. The stars are moving counterclockwise in the xy plane. 

Nonlinear Evolution of the Magnetothermal Instability in Two Dimensions by I. Parrish & J. Stone, ApJ 633, 334 (2005). Saturation of the Magnetothermal Instability in Three Dimensions by I. Parrish & J. Stone, ApJ 664, 135 (2007). Magnetic energy density along the faces of the computational volume in a threedimensional simulation of the MTI in a planeparallel stratified atmosphere. Magnetic field amplification by vertical convective motion is clearly evident. The MTI is driven by anisotropic conduction along magnetic field lines in a weakly collisional plasma, destabilizing atmospheres with entropy profiles that are stablystratified in the abscence of conduction. 

Nonlinear Evolution of the Magnetohydrodynamic RayleighTaylor Instability by J. Stone & T. Gardiner, Phys. Fluids, 19, 094104 (2007). The Magnetic RayleighTaylor Instability in Three Dimensions by J. Stone & T. Gardiner, ApJ 671, 1726 (2007).
Isosurfaces of the the density showing surfaces of the heavy (red) and
light (blue) fluids, as well as color slices of the density along the
edges of the computational domain, at late time in the hydrodynamic
RT instability (top) and strongly magnetized RT instability (bottom).
In the MHD case, the magnetic field is initially parallel to the xaxis.
The MHD flow shows much less fine scale structure than the hydro run,
due to suppression of modes parallel to B. Nonetheless, the interface
is still strongly unstable. 

The Magnetohydrodynamics of ShockCloud Interaction in Three Dimensions by M.S. Shin, J. Stone, & G. Snyder, ApJ 680, 336 (2008) Volumetric renderings of the density (left column) and magnetic energy (right column) after a Mach 10 shock interacts with a spherical magnetized cloud. The top row shows the case in which the magnetic field is parallel to the direction of shock propagation, the middle row shows the case in which it is oblique (at 45 degrees), and the bottom row shows the case in which it is perpendicular. In all cases the ratio of the initial gas to magnetic pressure is ten. 

Density Probability Distribution Functions in Supersonic Hydrodynamic and MHD Turbulence by N. Lemaster & J. Stone, ApJ 682, L97 (2008). Dissipation and Heating in Supersonic Hydrodynamic and MHD Turbulence by N. Lemaster & J. Stone, ApJ 691, 1092 (2009). Logarithm of the density (colors) on three faces of the computational volume and representative magnetic field lines (arrows) from both strong field MHD (top) and hydrodynamic (bottom) supersonic turbulence computed at a resolution of 1024^3. The sonic Mach number is about 10 in both cases, whereas the Alfvenic Mach number is 0.5 for the strong field MHD case. There are significant differences in the statistics of the density between the cases, evident from the images. Click on the images to download higher resolution versions. 

Buoyant Bubbles in Intracluster Gas: Effects of Magnetic Fields and Anisotropic Viscosity by R. Dong & J. Stone, ApJ 704, 1309 (2009). Slices of the density in the xy plane (left half of each plot) and yz plane (right half of each plot) at t=1 (left column), t=4 (middle column) and t=8 (right column) during the buoyant rise of a low density bubble in a sphericallysymmetric atmosphere with an initially horizontal magnetic field, and anisotropic viscosity relevant to low density plasmas. The top row shows the evolution of run with a weak field, the middle an intermediate strength field, and the bottom row a strong field. A linear color scale from d=0 to d=1 is used. 

Sustained Magnetorotational Turbulence in Local Simulations of Stratified Disks with Zero Net Magnetic Flux by S. Davis, J. Stone, & M. Pessah, ApJ 713, 52 (2010). Isosurface (at d=0.75) and slices of the density at 250 orbits in a domain of size 4H X 4H X 4H during the evolution of the magnetorotational instability in a vertically stratified disk. On the left face of the domain a slice of the magnitude of the magnetic field is shown. 