# Implications of thermodynamic inconsistency in CFD calculations

During my PhD work, I had to use tabulated values of thermodynamic properties of gases in some Computational Fluid Dynamics (CFD in short) simulations.

My tables are discretized in temperature and pressure (which are the independent variables) as reported in the scheme below: So I had to reconstruct an approximation of the general quantity $$\phi$$ for arbitrary couples of $$p$$,$$T$$ so that the CFD solver could retrieve them during the calculation.

Doing so with a simple Bilinear approach, works quite well for single species gas, but, when it comes to using mixtures I've observed really poor convergence rates and non-physical solutions.

The local bilinear approach it appears to be thermodynamic inconsistent: Maxwell relationships are not respected when thermodynamic quantities are independently interpolated.

On the other hand, using a consistent interpolation approach, such as a Hermite polynomial basis, everything works well without particular issues.

I would like to know if anyone has ever experienced such kind of problems, and if you have some ideas on the role of Maxwell relationships in non-ideal and compressible fluid dynamics.

• What are you using these arbitrary couples for? If it's visualization, it shouldn't make a difference. If it's something with the non linear solver, then it could very well be a problem – EMP Nov 18 '19 at 18:13
• My interpolations are used by the Roe approximate Riemann solver in order to provide the gas physical properties at any temperature a pressure. I would like to find an explication to what could be the problem with bilinear interpolations. I believe that thermodynamic consistency is involved: when I use a consistent approach (interpolations with Hermite polynomials) everything works well. I would like to understand why exactly. – iterrate Nov 19 '19 at 9:26
• Are you limiting the bilinear interpolation, to ensure that you don't overshoot on either face? If so, are you using face-based or cell-based limiting? Did you couple the bilinear approach so that you maintain a full concentration across all species during the interpolation? – EMP Nov 19 '19 at 15:49
• In the context of your question I think it is better if I clarify the figure in the main thread. The "cell" displayed is a schematization of the thermodynamic table structure. It has nothing to do with the spatial discretization of the fluid domain in which continuity, momentum and energy conservation are solved. Basically I had a big file that says (for instance) T: 300 K, P: 1 Bar -> Density : 1 kg/m3; T: 310 K, P: 1 bar, Density 2 kg ... and so on. Interpolation is required where the CFD solver wants to know what is the density value at T 302.254 K and P 1.325 bar. – iterrate Nov 19 '19 at 17:32
• Interesting question. I have a couple of questions in return: (1) You say that "bilinear approach it appears to be thermodynamic inconsistent" - Have you checked this by explicitly computing the various derivatives? (2) You say that "using a consistent interpolation approach, such as a Hermite polynomial basis, everything works well without particular issues" - Have you checked directly that the Maxwell relationships are satisfied or are you inferring this? – Brian Zatapatique Nov 23 '19 at 16:12

So the first thing I'm noticing is that if you have many variables, but they depend on each other through equations of state and other equations they are therefore coupled. I think it's best to interpolate the minimum number of variables required (maybe temperature and pressure) and then use your additional equations of state to calculate the other variables. Otherwise, you can end up with some nonphysical states when you interpolate every state. This should help keep your Roe states more feasible. Think about the example where you have a table of the values for the equation $$y=x^2$$, if you want to know the value at x=3.5, you should evaluate $$y=x^2$$ at x=3.5, not interpolate between the values of y at 3 and 4.
Furthermore, as you point out, the derivatives aren't continuous, maybe expand the stencil, and keep the interpolation $$C^1$$, so that any derivatives you have aren't jumping around like crazy depending on small changes in T. If you are using a strong solver (like Newton-Krylov), it uses the jacobian, and the derivatives of this interpolation would be jumping all over the place, as would the residual calculation (this can also be an issue with limiting the gradients in CFD solvers, and shows itself in the simulation never converging).
• Beside the computational expenses required to calculate the exact fluid composition (by minimising the Gibbs free enthalpy) are not feasible for the constraints that I have. The only things I can rely on are these tables, which have been calculated using accurate molecular dynamics simulations. I have tried smoother interpolations as bicubic, and Inverse distance weighting to provide at least $C^1$ functions: In such conditions I have "convergence" in the sense of low residuals, but I observe weird phenomenas occurring in my solutions (heat flux from colder zones to warm ones). – iterrate Nov 20 '19 at 22:29