
[Sponsors] 
June 3, 2020, 03:50 
psi in the rhoThermo class

#1 
Senior Member
Ruiyan Chen
Join Date: Jul 2016
Location: Hangzhou, China
Posts: 162
Rep Power: 9 
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 

June 9, 2020, 22:58 

#2 
Senior Member
Ruiyan Chen
Join Date: Jul 2016
Location: Hangzhou, China
Posts: 162
Rep Power: 9 
Anyone? Really looking forward to hearing your input here!


June 10, 2020, 04:43 

#3 
Senior Member

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 ddtterms. The code for pEqn in heat transfer/bouyantSimpleFoam contains fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA) * rhoSimpleFoam now instantiates the lowerlevel 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/energyequation/ 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 nonpremixed 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 pressuredensity 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 nonadiabatic case (e.g. in heat release by combustion)? * what is the role of the Mach number and the thermo dynamics (rho vs. psi thermodynamics)? Can we construct a high Mach number test case in which the rhothermodynamics 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, socalled wavetransmissive 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 xdirection, 1m in ydirection and 0.1 m in zdirection. 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/computationalaeroacousticspart2 [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/bcadvect...vetransmissive [8]: HowTo Using WaveTransmissive Boundary Conditions in OpenFoam: https://openfoamwiki.net/index.php/H...dary_condition 

June 10, 2020, 23:04 

#4  
Senior Member
Ruiyan Chen
Join Date: Jul 2016
Location: Hangzhou, China
Posts: 162
Rep Power: 9 
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 psibased one and a rhobased one, do you think the relation between rho and p, i.e., rho = psi * p, is still used in the rhobased 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 rhobased 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 psibased thermo, can you elaborate on the statement Quote:
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 rhobased thermo? I would think that wave is a typical example of "fluids with low compressibility" and it should be used with rhobased thermo. Thank you for taking the time to read my comments. 

June 15, 2020, 02:18 

#6 
Senior Member

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


June 15, 2020, 04:22 

#7 
Senior Member
Ruiyan Chen
Join Date: Jul 2016
Location: Hangzhou, China
Posts: 162
Rep Power: 9 
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 

June 15, 2020, 05:54 

#8 
Senior Member

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.


June 15, 2020, 10:37 

#9 
Senior Member

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. 

July 26, 2020, 07:32 

#10 
New Member
wangkuiming
Join Date: Jul 2020
Posts: 2
Rep Power: 0 
I plan to use rhoCentralFoam to simulate the liquid weak compressibility and thermal expansion, the state equation of ρ=ρ0[1+β(PP0)α(TT0)]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


October 21, 2020, 05:09 

#11 
Senior Member
Ruiyan Chen
Join Date: Jul 2016
Location: Hangzhou, China
Posts: 162
Rep Power: 9 
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 builtin 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. 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
"unknown psi type" error in reactingFoam  Speno93  OpenFOAM Running, Solving & CFD  2  May 25, 2021 09:15 
The udf.h headers are unable to open in VISUAL STUDIO 13  sanjeetlimbu  Fluent UDF and Scheme Programming  4  May 2, 2016 05:38 
Different define of psi uesd in pEqn.H  zqlhzx  OpenFOAM Running, Solving & CFD  0  December 24, 2013 08:45 
Nested class and inheritance permissions  MMC15  OpenFOAM Programming & Development  0  December 20, 2013 10:16 
Errors running allwmake in OpenFOAM141dev with WM_COMPILE_OPTION%3ddebug  unoder  OpenFOAM Installation  11  January 30, 2008 20:30 