[cryabroad]

Greetings Foamers,

As we know, there are psiThermo and rhoThermo classes to select from in FOAM. In psiThermo rho = psi*p. psi is given by the EOS, and by psi*p we get the density. My question is, in the rhoThermo class, is rho = psi*p still assumed? If yes, since both rho and psi are gived by the EOS and p is calculated, doesn't this mean that we overconstrained the system?

With this comes another more general question, why do we need both psiThermo and rhoThermo? I noticed that most heat transfer solvers use the rhoThermo class but most combustion ones use the psiThermo class, why the difference? I have a feeling that this might have something to do with the low Mach number assumption but couldn't quite figure it out.

Thanks in advance!
Ruiyan

[cryabroad]

Anyone? Really looking forward to hearing your input here!

[dlahaye]

Greetings,

1/

Thermodynamics in function of rho (density) or in terms of psi (compressibility)

* thermo.rho() returns a stored rho field that is calculated from pressure and temperature fields according to the selected thermophysical model in thermo.correct();
* rhothermo:*basic thermodynamic properties based on density; for fluids with low compressibility like water where the density change is mostly due to temperature variation; huge pressure jumps will not be represented correctly; solvers using rhothermo (like the buoyant solvers) use a simplified pressure equation where psi is accounted for explicitly;

Code for pEqn (as used for instance in heattransfer/buoyantPimpleFoam)
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA)
- fvm::laplacian(rAUf, p_rgh)
==
fvOptions(psi, p_rgh, rho.name())


* psiThermo:*basic thermodynamic properties based on compressibility; for gas; combustion solvers and transsonic solvers use the psithermo model. They have another pressure equation where the change in time of pressure is treated implicit by ddt(psi, p); larger pressure jumps are now possible. Such solvers are used when the flow is mainly driven by pressure changes and temperature change is only an effect of large compression and expansion;

Code for pEqn (as used for instance in combustion/reactingFoam)
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())

* note that alternatives above are essentially different in the ddt-terms. The code for pEqn in heat transfer/bouyantSimpleFoam contains
fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)

* rhoSimpleFoam now instantiates the lower-level fluidThermo which instantiates either a psiThermo or rhoThermo according to the 'type' specification in thermophysicalProperties. For details, see https://develop.openfoam.com/Develop...230a8965bdb860
* References: OpenFOAM.com: OpenFOAM: https://www.openfoam.com/documentati...mophysical.php
* References: CFD Direct: OpenFOAM: https://cfd.direct/openfoam/energy-equation/

2/ Notes in progress follow

Stationary vs. Transient Compressible Flow through a Rectangular Channel: Pressure Waves and the Wave Transmissive Boundary Condition

1/ Introduction

The goal of this humble endeavor is to complement existing documentation of compressible flow computations using the SIMPLE and PIMPLE algorithm. We target the simulation of turbulent non-premixed combustion in configurations similar to the RANS combustion tutorial SandiaD_LTS. The boundary condition setting for the pressure at the outlet patch will be looked into in particular. We are motivated by the study of turbulent combustion is a long cylindrical furnaces in which the inflow of air and oxidizer causes the flow to become compressible in a region close to the inlets. Our aim here is to systematically present information that we found in various sources.

This post consists of a theoretical and a numerical part. In the theoretical part we recall the derivation of the wave equation starting from the conservation of mass, the conservation of momentum and the pressure-density relationship for compressible flow. The adiabatic case is considered. Viscous and convective terms are neglected. We wish to make the point that the equations solved by rhoPimpleFoam (and by extension reactingFoam) are a generalization of the scalar wave equation. This theoretical part is unrelated to a particular CFD package.

In the numerical part we perform simulations using rhoSimpleFoam and rhoPimpleFoam. The geometry is deliberately kept simple to a rectangular channel. We wish to pursue two arguments. The first argument is the differences in convergence between rhoSimpleFoam and rhoPimpleFoam for this channel flow. This argument is relevant in understanding the difference in convergence of steady state and transient simulations for turbulent combustion in OpenFoam. An example of the former is the heat transfer tutorial reverseBurner. The second argument is the role of the wave transmissive outlet boundary conditions for the pressure in the case of transient simulations. This argument is important in understand when and how these boundary conditions should be applied.

2/ Derivation of the Wave Equation from the Compressible Flow Equations

It is well known that the scalar wave equation can be derived from the conservation of mass and momentum for compressible flow equation and the constitutive equation linking pressure and density. For the adiabatic case (i.e. without considering the conservation of energy and p = C \rho), this derivation is given in [1]. A slightly easier derivation assuming potential flow is given in [2]. This derivation is interesting to see for two reasons. The first is that is shows that velocity potential, pressure and density variations alike satisfy the wave equation. The second is that a velocity for the propagation of the wave is derived. This value allows to verify numerical simulations. Other derivations are given in [3,4]. More advanced derivations that include viscous effects are included in e.g. aeroacoustic as shown in [5]. Such viscous effects are included in rhoSimpleFoam and rhoPimpleFoam. 



Open Questions:
* what happens in the viscous case in the bulk flow (change molecular viscosity in transportProperties) and in the boundary layer (compare laminar and turbulent)? How does turbulent viscosity affect the viscous terms ? We can perform numerical simulations using rhoSimpleFoam and rhoPimpleFoam to find out;
* why are convective terms in the derivation systematically neglected? We can read the lecture notes by Rienstra and Hirshberg to find out more;
* what happens in the non-adiabatic case (e.g. in heat release by combustion)?
* what is the role of the Mach number and the thermo dynamics (rho vs. psi thermo-dynamics)? Can we construct a high Mach number test case in which the rho-thermodynamics clearly fails.

3/ Wave Transmissive Boundary Conditions

Acoustic waves reflect from hard walls. After hitting the wall, the wave changes sign and propagates back into the domain. This effect is beautifully illustrated in [6]. A Dirichlet or fixedValue (in OpenFoam) boundary conditions cause the boundary to become a hard wall. Due care is thus required in imposing fixed value boundary conditions. To prevent p to reflect from the wall, so-called wave-transmissive boundary conditions can be imposed. These boundary conditions are described in [5,6]. This boundary condition is used with the Sandia_LTS combustion RAS tutorial (!)

Open Questions:
* more details on how boundary conditions should be applied are required;
* a tutorial case/literature in which replacing waveTransmissive boundary conditions by fixedValue results in significantly different solution would be valuable to see;

4/ Rectangular Channel Case Definition

To be able to illustrate particular features, we decided to construct a test case. The geometry of the channel is defined as a rectangle measuring 5m in x-direction, 1m in y-direction and 0.1 m in z-direction.
The mesh is generated using blockMesh with setting specified in the blockMeshDict.
Derive wave speed.

5/ Incompressible Computations

In the incompressible case, the steady and transient simulations converge to the same steady state simulation. Steady state is elliptic. Transient is parabolic that converges to the steady state.

6/ Stationary Computation using rhoSimpleFoam

In case that stationarity is assumed, wave equation reduces to the Poisson equation. In the numerical simulations, no pressure reflections are seen;

7/ Transient Computation using rhoPimpleFoam

pressure waves with different when reflected from the outlet. Dependence on inlet velocity and molecular viscosity.

8/ Discussion and Remaining Questions

Remaining question:
numerics - not OpenFoam related: how does PIMPLE compare to Riemann solver;

9/ Conclusions

References
[1]: derivation of the wave equation for the linearized pressure for the adiabatic case with schematic representation: https://en.wikipedia.org/wiki/Acoustic_wave_equation
[2]: short derivation of the wave equation for the velocity potential in the adiabatic (no energy equation) and potential flow case: https://caefn.com/caa/computational-aeroacoustics-part2
[3]: extensive derivation of the wave equation for the linearized pressure for the adiabatic case: https://en.wikipedia.org/wiki/Acoustic_theory
[4]: book Rienstra and Hirchberg: search term “heat”, Section 5.5 Thermo Acoustics
[5]: aeroacoustic: https://en.wikipedia.org/wiki/Aeroacoustics
[6]: wave equation: https://en.wikipedia.org/wiki/Wave_equation
[7]: advective and waveTransmissive boundary conditions in OpenFoam: https://caefn.com/openfoam/bc-advect...vetransmissive
[8]: HowTo Using WaveTransmissive Boundary Conditions in OpenFoam: https://openfoamwiki.net/index.php/H...dary_condition

[cryabroad]

Hi Lahaye,

Thank you for sharing all these materials and references, some of them I have seen before but some of them definitely not. Up to now I only get a chance to browse and check part of them, will read more carefully later!

To your 1st point about thermo class, which is split to a psi-based one and a rho-based one, do you think the relation between rho and p, i.e., rho = psi * p, is still used in the rho-based thermo class? To me it is, as one can see from the psi*correction(fvm::ddt(p_rgh)) term in code. Following this, is it really rho = psi * p or rho = psi * p_rgh? More generally, I'm not sure why we still need psi in the rho-based thermo class. This psi*correction(fvm::ddt(p_rgh)) term should go to 0 at convergence, so it is added mainly due to numerical reasons I guess?

About the psi-based thermo, can you elaborate on the statement

They have another pressure equation where the change in time of pressure is treated implicit by ddt(psi, p); larger pressure jumps are now possible.

. I understand that by using rho = psi * p, ddt(rho) now becomes ddt(psi,p), which indicates implicit treatment of pressure, but why implicit treatments make larger pressure jumps possible? From a numerical point of view, it seems to me that the real reason we have rho-based and psi-based thermo is because of their different treatment of pressure. Basically, rhoThermo allows minor pressure jumps but psiThermo allows larger jumps?

To your 2nd point about waves, I'm actually very interested in such things. I see that both rhoSimpleFoam and rhoPimpleFoam are tested, what about solvers using rho-based thermo? I would think that wave is a typical example of "fluids with low compressibility" and it should be used with rho-based thermo.

Thank you for taking the time to read my comments.

[dlahaye]

Thank you for your reply. I need to take a closer look myself. Let us stay in touch. Best wishes.

[dlahaye]

I found this to be valuable: https://bugs.openfoam.org/view.php?id=1772#c5033

[cryabroad]

Hi Lahaye,

Thank you, I've read that before actually, it does provide very valuable information. I'm still looking for explanations to the idea that "implicit treatment of pressure allows for larger pressure jumps but explicit one doesn't", could you explain it in more detail or maybe pointing out some references?

Thanks,
Ruiyan

[dlahaye]

Not sure. I read however that rhoThermo is meant to handle low Mach numbers. I therefore imagine rhoThermo to break in case of high Mach number. One could e.g. take a rhoSimpleFoam tutorial case and see what happens in case that one gradually increases the mass flow rate at the inlet. Possibly one can place various simulations in a loop using foamDictionary.

[dlahaye]

rhoSimpleFoam is not limited to rhoThermo as explained here: https://develop.openfoam.com/Develop...230a8965bdb860

Not sure what is means to "break" rhoSimpleFoam using rhoTherm as imposing limits on either pressure or density or switching to rhoPimpleFoam might "fix" rhoThermo again.

Requires further thinking.

[wkm19920807]

I plan to use rhoCentralFoam to simulate the liquid weak compressibility and thermal expansion, the state equation of ρ=ρ0[1+β(P-P0)-α(T-T0)]to use include compression coefficient and coefficient of expansion, is constant, the expression with Boussinesq assumptions like but rhoCentralFoam only allow psiThermo type, more specifically the ideal gas (called perfectGas) and Peng Robinson EOS 1. can this change of state equation be realized, what documents are need to change?2. I know that Boussinesq assumes that it is in rhoThermo. What is the essential difference between psiThermo and rhoThermo?In psiThermo, I know that FOAM assumes rho = psi * p, and that psi is an ideal gas with a compressibility of 1/RT. If it were my equation of state in this form, what would the psi be

[cryabroad]

I came back with things I've learnt in the past several months and hopefully it sheds some lights on the original question I asked.

In the psiThermo class, psi = rho/p. The implication is that the EOS can be written in such a way that the density equals to something times the pressure. This, from the built-in EOS in FOAM, only applies to ideal gas and peng robinson gas. One can check this by using the rhoPimpleFoam solver (which uses psiThermo) and the "banana" method in the equationOfState entry in thermophysicalProperties file. This constraint has been relaxed in higher version FOAM though, because rhoPimpleFoam now allows the use of rhoThermo too.

In the rhoThermo class, psi should be the partial derivative of density with respect to pressure at constant temperature. This is clearly stated in Chapter 10.2 of Ferziger and Peric's book Computational Methods for Fluid Dynamics, 3rd ed.

The remaining question is still there, which is, why using two different thermo classes? When does one outperform the other? It seems to me that the rhoThermo class is a more general one and should be used for any solver.