# 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, 2019 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. Nov 19, 2019 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, 2019 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. Nov 19, 2019 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? Nov 23, 2019 at 16:12

We've been going back and forth in the comments and thanks for adding so much more information.

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

• First of all I want to thank you for the answer. I will give you some further explications on my problem. Unfortunately, I do not dispose of conventional equations of state in my specific case. The gas I'm studying undergoes to temperature and pressure levels for which molecular recombinations takes place. To give an example, CO2 at temperatures beyond 5000 K dissociates freeing CO, O2, C, O and ions. In such cases, EoS are not available. Nov 20, 2019 at 22:28
• 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). Nov 20, 2019 at 22:29
• Well as I said above, if you are interpolating every single quantity you use, you're almost certainly violating the relations between them. You're saying you can't calculate the relations, well then that's going to be a problem, because your interpolated states are physically feasible. In the jacobian matrix, assuming you've linearized the interpolation, you're going to see some weird jumps and discontinuities that are going to hinder the efficacy of the newton solver. And yea, the simpler the problem, the less of an issue you'll have with jumping and nonphysical states.
– EMP
Nov 21, 2019 at 15:12
• That should say that the interpolated states aren't physically feasible
– EMP
Nov 21, 2019 at 17:43
• This is a complicated question, but your solver can't differentiate between a physical and non-physical state. And how it affects your solver depends on how you coded it, so I can't answer that. But it can definitely be a problem. In limiting some people do extremum checks, and if the violation is small compared to the quantities its allowed to be unlimited. Maybe you can adapt that as you suggested for the Maxwell Relations, but you'll need to figure out a way to quantify that yourself, it's beyond my expertise.
– EMP
Nov 21, 2019 at 18:52