Ok, then I need to check if the definition of psi is correct in my solver. I'm using tabulated material models for everything. Right now, I have precalculated speed of sound to avoid calculating the derivative in each step, which is a bit expensive to calculate as it requires more table lookups. This solution is correct for my case, but if it's possible to keep the old formulation (with the same performance) it would be more elegant.

Yes and no :)

Yes, because you must implement psi as a coefficient between rho and p, rho = psi*p. In such assumption, terms d (rho) / dt `=` d(psi*p) / dt and d( rho U) / dx = d ( psi*U*p) / dx can be approximated with pisoCentralFoam scheme without changing pressure equation.

No, because in your formulation psi is function of pressure and this adds additional non-linearity to system. For example in my solver psi = 1 / (R*T) and it depends only on T.

Oh I just noticed that I made an error in my post, of course I implemented psi as rho(p, T) / (p + SMALL) and not the inverse.

I'm not sure I understood your post correctly. I agree that this kind of EOS introduces nonlinearity. The derivatives you talked about should still be valid with this EOS, as they don't derive in regards to p or rho.

Is this nonlinearity in determining the pressure already being solved by using multiple corrector iterations, or are additional terms necessary like in sonicLiquidFoam?

I guess I really need to spend some time trying to understand the derivation of the pressure, unfortunately I'm a one man show here :(

It means only, that it would be better if psi is not a function of pressure, f.e. for one component gas psi = f(T) = 1 / (R*T), for multicomponent gas psi = f(T, Yi), Yi - mass fractions and so on

Yes that's clear. But do you think I can get a correct solution with the current implementation? And would I need more corrector iterations compared to a linear EOS?

Today I have done some reading on the wiki and now I believe I understand the incompressible PISO algorithm and the derivation of the pressure equation, but I have not transferred the method to the compressible case yet to compare it with the source code.

I think, only practice can show whether this approach is applicable for your case

to get equation for pressure, do next things:
1) put EoS in continuity equations:

d(rho)/dt + d(rho U )/dx = d(psi*p)/dt + d(psi*p*U)/dx =
d(psi*p)/dt + d(psi*U*p)/dx = S_rho

2) relate change in velocity with pressure gradient and other terms:
A*U = H - d(p)/dx

3) put expression for velocity into continuity equation, formulated for pressure
d(psi*p)/dt + d( (psi*H/A) p )/dx - d ( (psi*p/A) * d(p)/dx )/dx = S_rho

Fluxes in the last equation multiplied by pressure are mass fluxes, satisfying continuity

Thanks for the hints, I will try to follow them.

In the meantime, I have done some tests on a real world application using my tabulated data with some inconclusive results. Please see the attached image. I have calculated an expanding flow using an initial high temperature, high pressure region at the top left of the geometry and calculated only flow and radiation. When I use 3 correctors I notice a change in the density, even in the region the shockwave hasn't reached yet. Any ideas about that?

It's somewhat difficult to tell which result is best. 2nd order Crank-Nicolson scheme didn't seem to have much of an influence (I went completely 2nd order). However, both changing the maxCo value and changing the number of correctors has an influence on the results. Computing time is somewhat limited for me so I would prefer to go with the maxCo=0.25 cases, but it seems that the time step error is still too large?

Edit: Forgot to mention that the mesh is completely orthogonal.

Your explanation together with this old thread helped me to understand the derivation of the pressure equation in compressible flow solvers like sonicFoam. This actually made me realize that I need the density source terms in the pressure equation, which I didn't put there yet :)

I think this also made me understand the problems involved with the nonlinearity in the material models here.
Basically, the pressure equation is not completely implicit anymore because the psi(p,T) dependence is handled explicitly, while the terms in pEqn are implicit. I still haven't completely understood the significance of this though. Does it mean that my time step size is limited by this explicit treatment? Or would multiple correctors be able to fix this nonlinearity? I still need to check at which point psi is reevaluated.

Ok I understood the problem some more now. Psi gets updated inside of thermo.correct(), which is called in hEqn.H in your solver, outside of the PISO loop. Inside the PISO loop, rho is calculated with rho=psi*p, using the new p, but the compressibility from the old pressure value of the previous time step. I think this is what is causing the wrong pressures in the results I showed above. This means that thermo.correct needs to be called inside of the PISO loop to update the compressibility after the pressure equation has been solved. It's still not clear to me if this is a solution that will converge properly, and I'm also not sure if I need to call it inside the nonorthogonal corrector loop (would have to update a few terms defined outside of this loop then, possibly hurting performance). thermo.correct is also comparitively expensive to call because more than the compressibility is updated in there, so it might be better to update only the compressibility somehow. This would require a different structure of the thermophysical libraries though.

I'm going to run a few calculations with this change and I'll report back if this fixes the density errors.

Hi, give me 2-3 days to write answer

Sure, you don't have to hurry :)
I have run some calculations in which I call thermo.correct() inside the PISO corrector loop and it seems to give consistent results when I use atleast 3 correctors. Using only 2 correctors gave me some differences in results. Everything above 2 correctors looks very similar. I think this can be regarded as a general flaw in many OpenFOAM solvers, as the thermophysical APIs explicitly allow a pressure dependence in the compressibility but most solvers won't be able to process this dependency properly. Do you think I should file a bug report?

Hi, i understand your question. In gas dynamics, we have at least 4 equation to solve: mass, momentum, energy and equation of state -

D(rho)/Dt = S_rho
D(rhoU)/Dt = S_U - dp/dx
D(rho*e)/Dt = S_e - d(p*U)/dx
p = rho*R*T

Implicit solvers of OpenFOAM are uses PISO method to solve this equations sequentially, for example:

after correcting velocity and pressure, you can evaluate addition terms for energy equation, like div(p*U) for internal energy or dp/dt for enthalpy, D(0.5*U*U)/Dt and others

If you want to couple energy with pressure in PISO loop, then, maybe you have to couple density and other equations (such as mass fractions) with new fluxes.

From my experience, it is better to couple only pressure and velocity with inner PISO iterations, while other variables can be coupled with outer iterations (PIMPLE).
In other words, all variables, that affects compressibility (psi) must be computed once before PISO loop.

That's exactly the point, the compressibility is affected by the pressure in my model, which gets modified inside the PISO loop. This means that compressibility needs to be calculated in each PISO loop iteration. Otherwise mass won't be conserved because rho is calculated at the end of the PISO loop using the compressibility. If this compressibility belongs to a wrong pressure then rho cannot be valid anymore.

Your solver (and other PISO based ones) currently do the following:
Instead the following is needed if psi=psi(p,T):

In this case, you can update only compressibility (psi) after each pressure equation step, without solving for new temperature. And you must check that your iterations are converging

That's basically what I would like to do, however this would require modifications in the thermophysical library, or a manual implementation of the compressibility model in the solver. To keep things simple I call thermo.correct() right now. This has some overhead because data table lookups are somewhat expensive, but with 3 corrector iterations per time step it's less than a second per time step.

Do you think I should file a bug report for other solvers regarding this issue? The implementation of the thermophysical library API in which psi has p and T parameters hints that the solvers should support such a dependency, while in reality they don't.

Well indeed if your PSI depends on p it should be better to correct the psi term after updating the pressure in each corrector. The problems are:

1) It can reduce the convergence rate if there are huge variations of psi with p

2) It increases the computational cost

Either way you can maintain 2 nCorrectors in the PISO reducing the timeStep

Fortunately my psi only depends slightly on p, so there is not much of a problem with convergence. I would rather stay at large time steps to get shorter calculation times. Using maxCo=0.25 with 3 correctors gives me results which only show a small temporal error. Depending on the case I have Mach numbers ranging between 0 and 3. Using this solver I get a speedup between 2-4x compared to a sonicFoam-based solver. I have also run some shockTube comparisons between this solver and sonicFoam and I get somewhat better results with this solver, so overall it seems to work better :)

I'm very pleased with your feed back. Can you tell about mesh quality, which you are using in computations?

Right now I'm using a simple mesh with orthogonal cells only. I have used surface-snapped meshes before but I used to get slower calculation speeds because the variance in cell sizes leads to higher maximum Courant numbers.
The mesh I'm using right now has about 400k cells.

The speedup I'm observing is mostly because this solver avoids the additional iterations per time step. I'm also calculating radiation using the P1 model, which is taking 50-60% of the calculation time per time step due to bad convergence. Avoiding the repeated calculation of the P1 equations helps the performance significantly here.