View Issue Details

IDProjectCategoryView StatusLast Update
0003740OpenFOAMBugpublic2021-10-18 16:03
Reporterpeksa Assigned To 
PrioritynormalSeveritymajorReproducibilityalways
Status newResolutionopen 
PlatformGNU/LinuxOSUbuntuOS Version15.04
Product Versiondev 
Summary0003740: Misleading pressure difference calculation in activePressureForceBaffleVelocityFvPatchVectorField.C
DescriptionDear developers,

In the boundary condition "activePressureForceBaffleVelocityFvPatchVectorField" a choice for opening a wall/cyclic baffle can be chosen to depend on either force or pressure difference.

However, the pressure difference is calculated as an explicit loop over patch faces with no weighting, leading to wrong behavior. See the problematic code snippet below:
-------------------------------------------------------------------------
            forAll(cyclicFaceCells, facei)
            {
                valueDiff += p[cyclicFaceCells[facei]];
            }
            forAll(nbrFaceCells, facei)
            {
                valueDiff -= p[nbrFaceCells[facei]];
            }
            Info<< "Pressure difference = " << valueDiff << endl;
        }
        if ((mag(valueDiff) > mag(minThresholdValue_)) || baffleActivated_)
        {
             ....
        }
-------------------------------------------------------------------------

One should use the same valueDiff definition as for the force but then divide the summed scalar by area to get a weighted difference:

            valueDiff = valueDiff/gSum(patch().magSf());


An example fix is presented earlier on ESI version:
https://develop.openfoam.com/Development/openfoam/-/commit/da6675803b167e7c51c1d2da4d6a0d95e297e1c1

Steps To ReproduceSee source code.
TagsNo tags attached.

Activities

peksa

2021-10-18 16:03

reporter   ~0012241

Something like this should work and showed sane performance in my initial simple test case.
activePressureForceBaffleVelocityFvPatchVectorField.patch (2,444 bytes)   
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
index 7060887..a4932a3 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C
@@ -240,37 +240,33 @@ void Foam::activePressureForceBaffleVelocityFvPatchVectorField::updateCoeffs()
         const labelList& nbrFaceCells = nbrPatch.patch().faceCells();
 
         scalar valueDiff = 0;
+        scalar patchArea = 0.0;
 
-        if (fBased_)
+         // Add this side
+        forAll(cyclicFaceCells, facei)
         {
-             // Add this side
-            forAll(cyclicFaceCells, facei)
-            {
-                valueDiff +=p[cyclicFaceCells[facei]]*mag(initCyclicSf_[facei]);
-            }
-
-            // Remove other side
-            forAll(nbrFaceCells, facei)
-            {
-                valueDiff -=p[nbrFaceCells[facei]]*mag(initCyclicSf_[facei]);
-            }
+            valueDiff +=p[cyclicFaceCells[facei]]*mag(initCyclicSf_[facei]);
+            patchArea += mag(initCyclicSf_[facei]);
+        }
 
+        // Remove other side
+        forAll(nbrFaceCells, facei)
+        {
+            valueDiff -=p[nbrFaceCells[facei]]*mag(initCyclicSf_[facei]);
+        }
+        
+        if(fBased_)
+        {
             Info<< "Force difference = " << valueDiff << endl;
         }
-        else // pressure based
+        else //pressure based then weighted by area
         {
-            forAll(cyclicFaceCells, facei)
-            {
-                valueDiff += p[cyclicFaceCells[facei]];
-            }
+            valueDiff = valueDiff/gSum(patch().magSf());
+            Info<< "Average pressure difference = " << valueDiff << endl;
+        }
 
-            forAll(nbrFaceCells, facei)
-            {
-                valueDiff -= p[nbrFaceCells[facei]];
-            }
+        reduce(valueDiff, sumOp<scalar>());
 
-            Info<< "Pressure difference = " << valueDiff << endl;
-        }
 
         if ((mag(valueDiff) > mag(minThresholdValue_)) || baffleActivated_)
         {

Issue History

Date Modified Username Field Change
2021-10-18 14:43 peksa New Issue
2021-10-18 16:03 peksa File Added: activePressureForceBaffleVelocityFvPatchVectorField.patch
2021-10-18 16:03 peksa Note Added: 0012241