Computational fluid dynamics (CFD) simulations involving non-ideal gases or gas mixtures rely on thermodynamic models which can be computationally expensive. Look-Up Tables (LUTs) have been extensively used to reduce such computational costs. During the simulation, the CFD solver accesses pre-computed discrete values of thermodynamic properties which are stored in a LUT. Interpolation techniques are then essential to accurately estimate the fluid local properties for any intermediate state between those available in the LUT. The interpolation technique employed needs to guarantee thermodynamic consistency, meaning the interpolated properties satisfy the fundamental laws of thermodynamics. To reconstruct an interpolation satisfying these constraints, bi-quintic Hermite polynomials have been proposed by Swesty [1], which are built from a thermodynamic potential and its successive derivatives. The same technique was successfully used by Zeman et al. [2] for laser plasma investigation. The present paper successfully applies the bi-quintic interpolation in the context of NICFD simulations and establishes its superior accuracy with respect to conventional bilinear interpolation.