# Simple finite volume method for compressible Navier-Stokes equations

I am interested in writing a simple, cell-centered, 2D FVM code for the unsteady, compressible Navier-Stokes equations (including shocks). Most of my experience is with finite difference and finite element methods. I did develop a finite volume code for Sod's problem as a learning exercise a while back. I've been searching through books and the literature for a simple algorithm which I can implement as a starting point, without much luck. Most of the more recent papers on finite volume methods are focused on higher-order accurate schemes, etc. I only really want second-order accuracy in space. Does anyone know of a good description of an algorithm for the compressible Navier-Stokes equations?

I'd like to avoid DG for now as shock capturing is still an active research area. Any insights or thoughts would be greatly appreciated.

• OpenFOAM uses a hybrid piso-simple scheme (referred to as a "pimple" scheme) to decouple the equations so that the can be solved and/or updated one at a time. A second order upwind scheme, combined with other 2nd order operator discretizations, should lead you to 2nd order accuracy overall. – Paul Oct 5 '16 at 14:19
• @Paul, the OpenFOAM documentation says that the "pimple" scheme is for incompressible flows. I'm interested in compressible, supersonic flows. – James Oct 5 '16 at 15:14
• Compressible solvers in OpenFoam are also built on the same principle of decoupling the equations and solving each one individually and looping over the equations until a convergence criterion is met. For example, see rhoPimpleFoam and sonicFoam. – Paul Oct 5 '16 at 15:46
• In my opinion the most simple would be to use $dU/dt=F(U)_x$ with $U=[\rho, \rho u, \rho E]$ and an equation of state for $p$. Then use explicit time integration with a second-order Runge-Kutta method and avoid the whole rho pimple stuff. – chris May 12 '17 at 14:40
• could you please contact me i need some help to formulate code of SIMPLE algorithm. zemamou.mounia@gmail.com – Mounia Z-m May 14 '18 at 21:08

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).