11

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:

enter image description here

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.

Original post in Physics SE.

Anton Menshov
  • 8,672
  • 7
  • 38
  • 94
iterrate
  • 111
  • 4
  • 2
    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 09:26
  • 1
    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
  • With such kind of structure, for a given couple $T$-$p$, it is possible to identify $T_i$, $T_{i+1}$, $p_i$ and $p_{i+1}$ such that:

    \begin{array}{l} T_i = \sup \left{ T_k \in T_g : T_k \leq T \right} \ \ T_{i+1} = \inf \left{ T_k \in T_g : T_k > T \right} \ \ p_i = \sup \left{ p_k \in p_g : p_k \leq p \right} \ \ p_{i+1} = \inf \left{ p_k \in p_g : p_k > p \right} \ \end{array}

    Where $T_g$ and $p_g$ are sets containing the temperature and pressure values for which the generic physical property $\phi$ is tabulated.

     – iterrate Nov 19 '19 at 17:44
  • Therefore, once I identify the 4 values $\phi_{i,j}$ (the nodes of the cell) I reconstruct a surface that passes for each of these points in a {\it local sense}. Which means that from one cell to another the representation is not the same. Using a bilinear approach, of course ensure the patched surfaces to be $C^0$, on the other hand derivatives are not continuous. I think that this may cause a bad precondinitioning of the matrix that has to be digested by the CFD solver. However, this is of course an hypothesis. – iterrate Nov 19 '19 at 17:55
  • Do you use the bilinear interpolation to approximate some derivatives, e.g., to compute sound speed ? – cfdlab Nov 20 '19 at 02:51
  • Exactly I interpolate density, enthalpy, entropy, specific heat at constant pressure, speed of sound, viscosity and thermal conductivity. – iterrate Nov 20 '19 at 08:18
  • 2
    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
  • Thank you Zapatique for your questions. It is really simple to verify the inconsistency of the bilinear approach. Without going in deep I'll give a simple example: – iterrate Nov 25 '19 at 08:18
  • $c_p=\frac{\partial h}{dT}p$. From bilinear interpolation one constructs: $\bar{h} = a_h + b_h T + c_h p + d_h Tp$ and $\bar{c_p} = a{c_p} + b_{c_p} T + c_{c_p} p + d_{c_p} Tp$. Deriving this representation of enthalpy with respect to temperature do not leads to the expression of $\bar{c_p}$. Hermite polynomials do that. I don't write down the basis here: the expressions are really long and tedious. I've performed CFD calculations with both methods, to verify what I'm saying. – iterrate Nov 25 '19 at 08:34

1 Answers1

2

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

EMP
  • 2,079
  • 10
  • 19
  • 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. – iterrate Nov 20 '19 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). – iterrate Nov 20 '19 at 22:29
  • This led me to the question about the thermodynamic stability subject: indeed one can find in litterature that violation of the second principle may occur if a thermodynamic potential do not respect specific convexity properties. This led me to use Hermite polynomials to ensure this properties. Indeed, everything works well with this kind of approach. But I'm still wondering what numerically happens in both case, in the Jacobian matrix. Nonetheless, bilinear interpolation still do the job if I use a single species gas. – iterrate Nov 20 '19 at 22:29
  • 1
    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 '19 at 15:12
  • 1
    That should say that the interpolated states aren't physically feasible – EMP Nov 21 '19 at 17:43
  • I totally agree with your statements. I may lack in numerical culture, but how do the solver "recognize" non-physical states and how it deals with them ?. Is there a metric for which one can state "ok, in this segment the interpolated values are violating Maxwell relations, but the discrepancies are small enough for the solution to show something meaningful" as I observe in single species calculation. Papers on this subject generally states that tables should be refined to reduce errors. – iterrate Nov 21 '19 at 18:00
  • 1
    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 '19 at 18:52
  • Thank you for the really delightful discussion about this fairly complicated subject :) – iterrate Nov 21 '19 at 19:06
  • Happy to try to help! – EMP Nov 26 '19 at 15:58