The numerical simulation of plasma flows including non-neutral regions is a considerable challenge due to the large variation in physical time scales involved: while electrons travel at velocities as high as millions of meters per second within the plasma sheaths and discharges, neutrals and ions travel at velocities 4-6 orders of magnitude lower. This large discrepancy between the time scales leads to what is commonly referred to as a "stiff" system of equations which require an excessive number of iterations to reach convergence in relation to the smoothness of the solution. We discuss recent advances in numerical methods that get rid of the stiffness within plasma systems of equations and enable the addition of detailed plasma sheath and discharge models to a Computational Fluid Dynamics (CFD) flow solver while not decreasing its computational efficiency.
