View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0004254OpenFOAMPatchpublic2025-07-09 13:262025-07-09 14:30
ReporterLiJixiao Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status closedResolutionsuspended 
PlatformGNU/LinuxOSOtherOS Version(please specify)
Product Version13 
Summary0004254: Small pressure difference is flooded by rounding error, consider gauge pressure
DescriptionFor compressible flows, it's common to represent pressure in absolute value. However when pressure drop is small, the pressure difference between cells is much less than absolute pressure. As a results, pressure difference is flooded by rounding error of floating point numbers. This makes pressure equation hard or unable to solve, and gives non-physical results. This can be fixed by allowing gauge pressure when solving p-U coupling and switch absolute pressure when updating physical properties.

A reproducible case is uploaded as `tube.tar.gz`. It simulates a heated laminar flow in square tube, the velocity is very small. Pressure drop is approximately 0.095Pa, much smaller than absolute value(1e5Pa). `fluid` is used as solver module. Both single precision and double precision gives non-physical pressure field, while single precision also has non-physical velocity field.

In chemical engineering process, absolute pressure can be terribly higher than 1atm, using gauge pressure will be helpful to those conditions.

I've made a patch to fix this issue. In this patch, `pOffset`was added to `rhoThermo` and `RhoFluidThermo` will compute physical properties with `p+pOffset`. `pOffset` will be optional for all thermophysical models based on `rhoThermo`. To use gauge pressure instead of absolute pressure, set `pOffset` and update fixed pressure defined in `0/p`. Pressure constraints should also be changed, if there is any. `pOffset` is optional and defaults to 0Pa to be compatible with absolute pressure cases.

Example:
```cpp
// File: 0/p
// Using gauge pressure
internalField uniform 0;
```
```cpp
// File: constant/physicalProperties
// Absolute pressure = p+pOffset
pOffset 1e5 [Pa];
```
Steps To Reproduce1. Decompress `tube.tar.gz`
2. run `blockMesh`
3. run `foamRun`
4. run `paraFoam`
TagsNo tags attached.

Activities

LiJixiao

2025-07-09 13:26

reporter  

0001-Add-pOffset-to-rhoThermo-to-support-gauge-pressure.patch (8,704 bytes)    
From 4257d65bf496c97c54f381de6c361b87d015282e Mon Sep 17 00:00:00 2001
From: ToKiNoBug <tokinobug@163.com>
Date: Wed, 9 Jul 2025 20:21:48 +0800
Subject: [PATCH] Add pOffset to rhoThermo to support gauge pressure

---
 .../basic/rhoFluidThermo/RhoFluidThermo.C     | 44 ++++++++++---------
 .../basic/rhoFluidThermo/rhoFluidThermo.H     |  1 +
 .../basic/rhoThermo/rhoThermo.C               | 20 ++++++++-
 .../basic/rhoThermo/rhoThermo.H               | 14 +++++-
 4 files changed, 55 insertions(+), 24 deletions(-)

diff --git a/src/thermophysicalModels/basic/rhoFluidThermo/RhoFluidThermo.C b/src/thermophysicalModels/basic/rhoFluidThermo/RhoFluidThermo.C
index fec017b057..f304b74586 100644
--- a/src/thermophysicalModels/basic/rhoFluidThermo/RhoFluidThermo.C
+++ b/src/thermophysicalModels/basic/rhoFluidThermo/RhoFluidThermo.C
@@ -32,6 +32,8 @@ void Foam::RhoFluidThermo<BaseThermo>::calculate()
 {
     const scalarField& hCells = this->he();
     const scalarField& pCells = this->p_;
+    // Add pressureOffset to p, allowing gauge pressure
+    const scalar pressureOffset = this->pOffset().value();
 
     scalarField& TCells = this->T_.primitiveFieldRef();
     scalarField& CpCells = this->Cp_.primitiveFieldRef();
@@ -57,18 +59,18 @@ void Foam::RhoFluidThermo<BaseThermo>::calculate()
         TCells[celli] = thermoMixture.The
         (
             hCells[celli],
-            pCells[celli],
+            pCells[celli]+pressureOffset,
             TCells[celli]
         );
 
-        CpCells[celli] = thermoMixture.Cp(pCells[celli], TCells[celli]);
-        CvCells[celli] = thermoMixture.Cv(pCells[celli], TCells[celli]);
-        psiCells[celli] = thermoMixture.psi(pCells[celli], TCells[celli]);
-        rhoCells[celli] = thermoMixture.rho(pCells[celli], TCells[celli]);
+        CpCells[celli] = thermoMixture.Cp(pCells[celli]+pressureOffset, TCells[celli]);
+        CvCells[celli] = thermoMixture.Cv(pCells[celli]+pressureOffset, TCells[celli]);
+        psiCells[celli] = thermoMixture.psi(pCells[celli]+pressureOffset, TCells[celli]);
+        rhoCells[celli] = thermoMixture.rho(pCells[celli]+pressureOffset, TCells[celli]);
 
-        muCells[celli] = transportMixture.mu(pCells[celli], TCells[celli]);
+        muCells[celli] = transportMixture.mu(pCells[celli]+pressureOffset, TCells[celli]);
         kappaCells[celli] =
-            transportMixture.kappa(pCells[celli], TCells[celli]);
+            transportMixture.kappa(pCells[celli]+pressureOffset, TCells[celli]);
     }
 
     volScalarField::Boundary& pBf =
@@ -124,15 +126,15 @@ void Foam::RhoFluidThermo<BaseThermo>::calculate()
                     transportMixture =
                     this->transportMixture(composition, thermoMixture);
 
-                phe[facei] = thermoMixture.he(pp[facei], pT[facei]);
+                phe[facei] = thermoMixture.he(pp[facei]+pressureOffset, pT[facei]);
 
-                pCp[facei] = thermoMixture.Cp(pp[facei], pT[facei]);
-                pCv[facei] = thermoMixture.Cv(pp[facei], pT[facei]);
-                ppsi[facei] = thermoMixture.psi(pp[facei], pT[facei]);
-                prho[facei] = thermoMixture.rho(pp[facei], pT[facei]);
+                pCp[facei] = thermoMixture.Cp(pp[facei]+pressureOffset, pT[facei]);
+                pCv[facei] = thermoMixture.Cv(pp[facei]+pressureOffset, pT[facei]);
+                ppsi[facei] = thermoMixture.psi(pp[facei]+pressureOffset, pT[facei]);
+                prho[facei] = thermoMixture.rho(pp[facei]+pressureOffset, pT[facei]);
 
-                pmu[facei] = transportMixture.mu(pp[facei], pT[facei]);
-                pkappa[facei] = transportMixture.kappa(pp[facei], pT[facei]);
+                pmu[facei] = transportMixture.mu(pp[facei]+pressureOffset, pT[facei]);
+                pkappa[facei] = transportMixture.kappa(pp[facei]+pressureOffset, pT[facei]);
             }
         }
         else
@@ -149,15 +151,15 @@ void Foam::RhoFluidThermo<BaseThermo>::calculate()
                     transportMixture =
                     this->transportMixture(composition, thermoMixture);
 
-                pT[facei] = thermoMixture.The(phe[facei], pp[facei], pT[facei]);
+                pT[facei] = thermoMixture.The(phe[facei], pp[facei]+pressureOffset, pT[facei]);
 
-                pCp[facei] = thermoMixture.Cp(pp[facei], pT[facei]);
-                pCv[facei] = thermoMixture.Cv(pp[facei], pT[facei]);
-                ppsi[facei] = thermoMixture.psi(pp[facei], pT[facei]);
-                prho[facei] = thermoMixture.rho(pp[facei], pT[facei]);
+                pCp[facei] = thermoMixture.Cp(pp[facei]+pressureOffset, pT[facei]);
+                pCv[facei] = thermoMixture.Cv(pp[facei]+pressureOffset, pT[facei]);
+                ppsi[facei] = thermoMixture.psi(pp[facei]+pressureOffset, pT[facei]);
+                prho[facei] = thermoMixture.rho(pp[facei]+pressureOffset, pT[facei]);
 
-                pmu[facei] = transportMixture.mu(pp[facei], pT[facei]);
-                pkappa[facei] = transportMixture.kappa(pp[facei], pT[facei]);
+                pmu[facei] = transportMixture.mu(pp[facei]+pressureOffset, pT[facei]);
+                pkappa[facei] = transportMixture.kappa(pp[facei]+pressureOffset, pT[facei]);
             }
         }
     }
diff --git a/src/thermophysicalModels/basic/rhoFluidThermo/rhoFluidThermo.H b/src/thermophysicalModels/basic/rhoFluidThermo/rhoFluidThermo.H
index b0da199e72..5f69c055a1 100644
--- a/src/thermophysicalModels/basic/rhoFluidThermo/rhoFluidThermo.H
+++ b/src/thermophysicalModels/basic/rhoFluidThermo/rhoFluidThermo.H
@@ -55,6 +55,7 @@ namespace Foam
 class rhoFluidThermo
 :
     virtual public rhoThermo,
+    //virtual public rhoThermo,
     virtual public fluidThermo
 {
 public:
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
index 2403b915ce..324589599c 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
@@ -18,7 +18,7 @@ License
     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
     for more details.
 
-    You should have received a copy of the GNU General Public License
+    You should have received a copy of the GNU General Public Licensedisclose
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 \*---------------------------------------------------------------------------*/
@@ -46,8 +46,17 @@ Foam::rhoThermo::implementation::implementation
         ),
         mesh,
         dimDensity
+    ),
+    pOffset_(
+        dict.lookupOrDefault<dimensionedScalar>
+        (
+            "pOffset",
+            dimensionedScalar{"pOffset",dimPressure,0.0}
+        )
     )
-{}
+{
+    Info<<"Pressure offset is "<<this->pOffset().value()<<" Pa\n";
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -83,5 +92,12 @@ Foam::volScalarField& Foam::rhoThermo::implementation::rho()
     return rho_;
 }
 
+const Foam::dimensionedScalar & Foam::rhoThermo::implementation::pOffset() const {
+    return this->pOffset_;
+}
+
+Foam::dimensionedScalar & Foam::rhoThermo::implementation::pOffset() {
+    return this->pOffset_;
+}
 
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
index c12d3831be..3193b04933 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
@@ -77,6 +77,12 @@ public:
 
             //- Return non-const access to the local density field [kg/m^3]
             virtual volScalarField& rho() = 0;
+
+            //- Pressure offset when correcting physical properties [Pa]
+            virtual const dimensionedScalar &pOffset() const = 0;
+
+            //- Return non-const access to pressure offset [Pa]
+            virtual dimensionedScalar& pOffset() = 0;
 };
 
 
@@ -94,7 +100,8 @@ protected:
 
         //- Density field [kg/m^3]
         volScalarField rho_;
-
+        //- Pressure offset when correcting physical properties [Pa]. 0 [Pa] by default.
+        dimensionedScalar pOffset_;
 
 public:
 
@@ -124,6 +131,11 @@ public:
             //- Return non-const access to the local density field [kg/m^3]
             virtual volScalarField& rho();
 
+            //- Pressure offset when correcting physical properties [Pa]
+            virtual const dimensionedScalar &pOffset() const;
+
+            //- Return non-const access to pressure offset [Pa]
+            virtual dimensionedScalar& pOffset();
 
     // Member Operators
 
-- 
2.50.0
0001-Add-pOffset-to-rhoThermo-to-support-gauge-pressure.patch (8,704 bytes)   
p.png (272,175 bytes)
U.png (288,150 bytes)
tube.tar.gz (1,898 bytes)

henry

2025-07-09 13:50

manager   ~0013591

https://openfoam.org/maintenance/
https://openfoam.org/news/funding-2025/

henry

2025-07-09 13:52

manager   ~0013592

We can add support for pressure referencing in all the compressible solvers but we would not do it in the manner you suggest.
We will require maintenance funding for this work.

LiJixiao

2025-07-09 14:07

reporter   ~0013593

Thank you for considering gauge pressure!

I understand that my patch is far from perfect. I'm willing to contribute according to your design. Look forward to your opinions on pressure referencing(gauge pressure).

henry

2025-07-09 14:11

manager   ~0013594

Will you be providing maintenance funding to cover our work on this and future maintenance work?

LiJixiao

2025-07-09 14:24

reporter   ~0013595

I'm sorry to say that even silver plan is far beyond my ability. I'm just a PhD student that you can see everywhere in every college. However I'm still wiling to contribute code according to your design.

henry

2025-07-09 14:30

manager   ~0013596

Pending maintenance funding

Issue History

Date Modified Username Field Change
2025-07-09 13:26 LiJixiao New Issue
2025-07-09 13:26 LiJixiao File Added: 0001-Add-pOffset-to-rhoThermo-to-support-gauge-pressure.patch
2025-07-09 13:26 LiJixiao File Added: p.png
2025-07-09 13:26 LiJixiao File Added: U.png
2025-07-09 13:26 LiJixiao File Added: tube.tar.gz
2025-07-09 13:50 henry Note Added: 0013591
2025-07-09 13:52 henry Note Added: 0013592
2025-07-09 14:07 LiJixiao Note Added: 0013593
2025-07-09 14:11 henry Note Added: 0013594
2025-07-09 14:24 LiJixiao Note Added: 0013595
2025-07-09 14:30 henry Assigned To => henry
2025-07-09 14:30 henry Status new => closed
2025-07-09 14:30 henry Resolution open => suspended
2025-07-09 14:30 henry Note Added: 0013596