There are many ways to go about solving the compressible navier stokes equations. One approach is to solve all of the equations (continuity, momentum, & energy) together as one system of coupled equations. If your flow is inviscid, this is a fairly simple system of equations to solve. But if you want to model viscous flow, your resulting system of non-linear equations will be extremely large and may not even fit into your machine's memory. Furthermore, solving a system of non-linear equations presents its own set of challenges unto itself.
An alternative approach to resolving viscous flow is to decouple the equations and solve each, in a one-at-a-time fashion over an iterative loop within each timestep. This allows to you work more flexibly with limited computer memory constraints. It also linearizes the advection term in the momentum equation (and the kinetic energy term in the energy equation) so that you don't need to solve a non-linear system for velocity or energy. The PISO (pressure implicit splitting of operators) and SIMPLE (Semi-Implicit Method for Pressure Linked Equations) are both based upon this principle. Though PISO & SIMPLE were originally formulated to solve the incompressible navier stokes equations, they can also be used to solve the compressible version as well.
A hybrid between these two algorithms is sometimes referred to by the as the "PIMPLE" scheme (a portmanteau of piso and simple). This hybrid algorithm works as follows:
for each timestep
for a fixed number of iterations N1
solve continuity equation (for density)
solve momentum predictor equation (for velocity)
solve energy equation (for energy)
for some fixed number of iterations N2
solve pressure equation
update velocity
update density (using the equation of state)
end
solve turbulence equations
end
end
At any one given step, you use the most recently resolved quantities. For example, after solving the continuity equation for density, you use this density solution to solve the momentum predictor equation for velocity, then use this velocity to solve the energy equation, etc...
You will notice that this algorithm has both an outer loop (inspired by the SIMPLE algorithm) and an inner loop (inspired by the PISO algorithm). This is how the compressible navier stokes equations are resolved in OpenFOAM's rhoPimpleFoam solver.
Accuracy depends on your mesh/timestep sizes and on how you discretize the operators in each equation. For time derivatives, you can use a crank nicholson scheme to acheive second order accuracy in time. Using centered difference discretization for laplacian operators (which appear in the momentum predictor, energy, and pressure equations) will result in second order accuracy. However, the continuity, momentum predictor, and energy equations all contain 1st order divergence operators (think of them like first derivatives in space). A centered difference on a divergence operator is unstable. You might be able to achieve second order accuracy by using a second order upwind scheme for your divergence operators. Depending on the complexity of your flow (e.g. shock waves and other high mach compressible flow effects), a second order scheme may not result in second order accuracy and/or may not even be stable. The safest thing to do in such cases is use a first order upwind scheme (even though it lowers your overall accuracy, it generally leads to more stable time-stepping and/or shock capturing).