ENH: added surfaceFieldValue uniformity operation

src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2017-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,18 +83,21 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::operationTypeNames_
     { operationType::opCoV, "CoV" },
     { operationType::opAreaNormalAverage, "areaNormalAverage" },
     { operationType::opAreaNormalIntegrate, "areaNormalIntegrate" },
+    { operationType::opUniformity, "uniformity" },
 
     // Using weighting
     { operationType::opWeightedSum, "weightedSum" },
     { operationType::opWeightedAverage, "weightedAverage" },
     { operationType::opWeightedAreaAverage, "weightedAreaAverage" },
     { operationType::opWeightedAreaIntegrate, "weightedAreaIntegrate" },
+    { operationType::opWeightedUniformity, "weightedUniformity" },
 
     // Using absolute weighting
     { operationType::opAbsWeightedSum, "absWeightedSum" },
     { operationType::opAbsWeightedAverage, "absWeightedAverage" },
     { operationType::opAbsWeightedAreaAverage, "absWeightedAreaAverage" },
     { operationType::opAbsWeightedAreaIntegrate, "absWeightedAreaIntegrate" },
+    { operationType::opAbsWeightedUniformity, "absWeightedUniformity" },
 };
 
 const Foam::Enum
@@ -117,10 +120,8 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::obr() const
     {
         return mesh_.lookupObject<objectRegistry>(regionName_);
     }
-    else
+
-    {
+    return mesh_;
-        return mesh_;
-    }
 }
 
 
@@ -599,7 +600,7 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::initialise
         }
     }
 
-    // Backwards compatibility for v1612+ and older
+    // Backwards compatibility for v1612 and older
     List<word> orientedFields;
     if (dict.readIfPresent("orientedFields", orientedFields))
     {
@@ -696,6 +697,45 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
 
             return gSum(pos0(nv)*mag(values) - neg(nv)*mag(values));
         }
+
+        case opUniformity:
+        case opWeightedUniformity:
+        case opAbsWeightedUniformity:
+        {
+            const scalar areaTotal = gSum(mag(Sf));
+            tmp<scalarField> areaVal(values * mag(Sf));
+
+            scalar mean, numer;
+
+            if (canWeight(weightField))
+            {
+                // Weighted quantity = (Weight * phi * dA)
+
+                tmp<scalarField> weight(weightingFactor(weightField));
+
+                // Mean weighted value (area-averaged)
+                mean = gSum(weight()*areaVal()) / areaTotal;
+
+                // Abs. deviation from weighted mean value
+                numer = gSum(mag(weight*areaVal - (mean * mag(Sf))));
+            }
+            else
+            {
+                // Unweighted quantity = (1 * phi * dA)
+
+                // Mean value (area-averaged)
+                mean = gSum(areaVal()) / areaTotal;
+
+                // Abs. deviation from mean value
+                numer = gSum(mag(areaVal - (mean * mag(Sf))));
+            }
+
+            // Uniformity index
+            const scalar ui = 1 - numer/(2*mag(mean*areaTotal) + ROOTVSMALL);
+
+            return min(max(ui, 0), 1);
+        }
+
         default:
         {
             // Fall through to other operations
@@ -742,6 +782,45 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::processValues
             const scalar val = gSum(values & Sf);
             return vector(val, 0, 0);
         }
+
+        case opUniformity:
+        case opWeightedUniformity:
+        case opAbsWeightedUniformity:
+        {
+            const scalar areaTotal = gSum(mag(Sf));
+            tmp<scalarField> areaVal(values & Sf);
+
+            scalar mean, numer;
+
+            if (canWeight(weightField))
+            {
+                // Weighted quantity = (Weight * phi . dA)
+
+                tmp<scalarField> weight(weightingFactor(weightField));
+
+                // Mean weighted value (area-averaged)
+                mean = gSum(weight()*areaVal()) / areaTotal;
+
+                // Abs. deviation from weighted mean value
+                numer = gSum(mag(weight*areaVal - (mean * mag(Sf))));
+            }
+            else
+            {
+                // Unweighted quantity = (1 * phi . dA)
+
+                // Mean value (area-averaged)
+                mean = gSum(areaVal()) / areaTotal;
+
+                // Abs. deviation from mean value
+                numer = gSum(mag(areaVal - (mean * mag(Sf))));
+            }
+
+            // Uniformity index
+            const scalar ui = 1 - numer/(2*mag(mean*areaTotal) + ROOTVSMALL);
+
+            return vector(min(max(ui, 0), 1), 0, 0);
+        }
+
         default:
         {
             // Fall through to other operations
@@ -762,11 +841,9 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor
     {
         return mag(weightField);
     }
-    else
+
-    {
+    // pass through
-        // pass through
+    return weightField;
-        return weightField;
-    }
 }
 
 
@@ -788,10 +865,8 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor
     {
         return mag(weightField * mag(Sf));
     }
-    else
+
-    {
+    return (weightField * mag(Sf));
-        return (weightField * mag(Sf));
-    }
 }
 
 
@@ -813,10 +888,8 @@ Foam::functionObjects::fieldValues::surfaceFieldValue::weightingFactor
     {
         return mag(weightField & Sf);
     }
-    else
+
-    {
+    return (weightField & Sf);
-        return (weightField & Sf);
-    }
 }
 
 

src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H

@@ -145,6 +145,21 @@ Note
         - take care when using isoSurfaces - these might have duplicate
           triangles and so integration might be wrong
 
+    Uniformity
+    \f[
+        UI(\phi) = 1 - \frac{1}{2 \overline{\phi} A}
+        \int{\left| W \phi \cdot \hat{n} - \bar{W} \bar{\phi}\right| d\vec{A}}
+        \,,\;
+        \bar{\phi} = \frac{\int{W \phi \cdot d\vec{A}}}{\int{W \cdot d\vec{A}}}
+    \f]
+
+    A velocity uniformity index is calculated with no weighting (W=1) and
+    \f$ \phi = \vec{U} \f$.
+
+    A scalar concentration uniformity index is calculated with either
+    \f$ \rho \vec{U} \f$ or \f$ \vec{U} \f$ for weighting and
+    \f$ \phi = conc \f$.
+
 See also
     Foam::fieldValues
     Foam::functionObject
@@ -232,6 +247,7 @@ public:
             opCoV,                  //!< Coefficient of variation
             opAreaNormalAverage,    //!< Area average in normal direction
             opAreaNormalIntegrate,  //!< Area integral in normal direction
+            opUniformity,           //!< Uniformity index
 
             // Weighted variants
 
@@ -247,6 +263,9 @@ public:
             //! Weighted area integral
             opWeightedAreaIntegrate = (opAreaIntegrate | typeWeighted),
 
+            //! Weighted uniformity index
+            opWeightedUniformity = (opUniformity | typeWeighted),
+
             // Variants using absolute weighting
 
             //! Sum using abs weighting
@@ -261,6 +280,10 @@ public:
             //! Area integral using abs weighting
             opAbsWeightedAreaIntegrate =
                 (opWeightedAreaIntegrate | typeAbsolute),
+
+            //! Uniformity index using abs weighting
+            opAbsWeightedUniformity =
+                (opWeightedUniformity | typeAbsolute),
         };
 
         //- Operation type names

src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C

@@ -275,10 +275,28 @@ processSameTypeValues
 
         case opAreaNormalAverage:
         case opAreaNormalIntegrate:
+        case opUniformity:
         {
             // Handled in specializations only
             break;
         }
+
+        case opWeightedUniformity:
+        case opAbsWeightedUniformity:
+        {
+            if (canWeight(weightField))
+            {
+                // Change weighting from vector -> scalar and dispatch again
+                return processValues<Type, scalar>
+                (
+                    values,
+                    Sf,
+                    weightingFactor(weightField)
+                );
+            }
+
+            break;
+        }
     }
 
     return result;

tutorials/compressible/rhoSimpleFoam/squareBend/0/T

@@ -20,23 +20,22 @@ internalField   uniform 1000;
 
 boundaryField
 {
-    Default_Boundary_Region
-    {
-        type            zeroGradient;
-    }
-
     inlet
     {
         type            fixedValue;
-        value           uniform 1000;
+        value           $internalField;
     }
 
     outlet
     {
         type            inletOutlet;
-        //type            zeroGradient;
+        value           $internalField;
-        value           uniform 1000;
+        inletValue      $internalField;;
-        inletValue      uniform 1000;
+    }
+
+    ".*"
+    {
+        type            zeroGradient;
     }
 }
 

tutorials/compressible/rhoSimpleFoam/squareBend/0/U

@@ -33,8 +33,8 @@ boundaryField
     outlet
     {
         type            inletOutlet;
-        value           uniform (0 0 0);
+        value           $internalField;
-        inletValue      uniform (0 0 0);
+        inletValue      $internalField;
     }
 }
 

tutorials/compressible/rhoSimpleFoam/squareBend/0/alphat

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      alphat;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -25,17 +24,17 @@ boundaryField
     {
         type            compressible::alphatWallFunction;
         Prt             0.85;
-        value           uniform 0;
+        value           $internalField;
     }
     inlet
     {
         type            calculated;
-        value           uniform 0;
+        value           $internalField;
     }
     outlet
     {
         type            calculated;
-        value           uniform 0;
+        value           $internalField;
     }
 }
 

tutorials/compressible/rhoSimpleFoam/squareBend/0/epsilon

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      epsilon;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -27,19 +26,19 @@ boundaryField
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
-        value           uniform 200;
+        value           $internalField;
     }
     inlet
     {
         type            turbulentMixingLengthDissipationRateInlet;
         mixingLength    0.005;
-        value           uniform 200;
+        value           $internalField;
     }
     outlet
     {
         type            inletOutlet;
-        inletValue      uniform 200;
+        value           $internalField;
-        value           uniform 200;
+        inletValue      $internalField;
     }
 }
 

tutorials/compressible/rhoSimpleFoam/squareBend/0/k

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      k;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -24,19 +23,19 @@ boundaryField
     Default_Boundary_Region
     {
         type            kqRWallFunction;
-        value           uniform 1;
+        value           $internalField;
     }
     inlet
     {
         type            turbulentIntensityKineticEnergyInlet;
         intensity       0.05;
-        value           uniform 1;
+        value           $internalField;
     }
     outlet
     {
         type            inletOutlet;
-        inletValue      uniform 1;
+        value           $internalField;
-        value           uniform 1;
+        inletValue      $internalField;
     }
 }
 

tutorials/compressible/rhoSimpleFoam/squareBend/0/nut

@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0";
     object      nut;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -27,17 +26,17 @@ boundaryField
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
-        value           uniform 0;
+        value           $internalField;
     }
     inlet
     {
         type            calculated;
-        value           uniform 0;
+        value           $internalField;
     }
     outlet
     {
         type            calculated;
-        value           uniform 0;
+        value           $internalField;
     }
 }
 

tutorials/compressible/rhoSimpleFoam/squareBend/0/p

@@ -28,17 +28,17 @@ boundaryField
     {
         type            zeroGradient;
         //type            mixed;
-        refValue        uniform 110000;
+        refValue        $internalField;
         refGradient     uniform 0;
         valueFraction   uniform 0.3;
     }
     outlet
     {
         type            fixedValue;
-        value           uniform 110000;
+        value           $internalField;
 
         //type            mixed;
-        //refValue        uniform 110000;
+        //refValue        $internalField;
         //refGradient     uniform 0;
         //valueFraction   uniform 1;
         //type            transonicOutletPressure;
@@ -46,7 +46,12 @@ boundaryField
         //phi             phi;
         //gamma           1.4;
         //psi             psi;
-        //pInf            uniform 110000;
+        //pInf            $internalField;
+    }
+
+    ".*"
+    {
+        type            zeroGradient;
     }
 }
 

tutorials/compressible/rhoSimpleFoam/squareBend/system/controlDict

@@ -47,5 +47,12 @@ graphFormat     raw;
 
 runTimeModifiable true;
 
+#include "sampleControls"
+
+functions
+{
+    #include "sampling"
+}
+
 
 // ************************************************************************* //

tutorials/compressible/rhoSimpleFoam/squareBend/system/fieldTransfer

@@ -0,0 +1,51 @@
+// -*- C++ -*-
+
+// ************************************************************************* //
+
+// Transcribe volume fields to surfaces.
+fieldTransfer
+{
+    type    surfMeshes;
+    libs    ("libsampling.so");
+    log     true;
+    writeControl    none;
+    createOnRead    true;
+    executeControl  timeStep;
+    executeInterval 1;
+
+    fields      (p rho U T);
+    derived     (rhoU);
+
+    _plane
+    {
+        type    plane;
+        source  cells;
+
+        planeType   pointAndNormal;
+
+        pointAndNormalDict
+        {
+            normal (-1 0 0);
+            point  (-0.04 0 0);
+        }
+    }
+
+    surfaces
+    (
+        // Top channel
+        plane1
+        {
+            ${_plane}
+            bounds (-1 0 -1) (0 1 1);
+        }
+
+        // Bottom channel
+        plane2
+        {
+            ${_plane}
+            bounds (-1 -1 -1) (0 0 1);
+        }
+    );
+}
+
+// ************************************************************************* //

tutorials/compressible/rhoSimpleFoam/squareBend/system/sampleControls

@@ -0,0 +1,40 @@
+// -*- C++ -*-
+
+//   restartTime:
+//     - a 'one-shot' reset at a particular time
+//
+// fields [required]
+//     Pairs of fields to use for calculating the deviation.
+//     The fields must already exist on the surfaces.
+//
+// weightField [optional]
+//     A scalar or vector field for weighting.
+//
+// postOperation [optional]
+//     Modify the results by particular operations.
+//     (none | sqrt)
+//     The sqrt operation is useful when determining RMS values.
+//
+// The 'output/write' control triggers the calculation.
+__surfaceFieldValue
+{
+    type            surfaceFieldValue;
+    libs            ("libfieldFunctionObjects.so");
+    log             on;
+    enabled         true;
+
+    writeControl    timeStep;
+    writeInterval   1;
+
+    writeFields     false;
+    surfaceFormat   vtk;
+    // writeArea       true;
+
+    // resetOnStartUp  true;
+    // resetOnOutput   false;
+    // periodicRestart true;
+    // restartPeriod   0.0005;
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling

@@ -0,0 +1,99 @@
+// -*- C++ -*-
+
+// ************************************************************************* //
+
+#include "fieldTransfer"
+
+massflow
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane1;
+
+    operation       areaNormalIntegrate;
+
+    fields          ( rhoU );
+}
+
+areaIntegrate
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane1;
+
+    operation       weightedAreaIntegrate;
+    weightField     rhoU;
+    fields          ( T );
+}
+
+// Inflow uniformity
+UI1
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane1;
+
+    operation       uniformity;
+    fields          ( U T );
+}
+
+
+// Uniformity after the bend
+UI2
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane2;
+
+    operation       uniformity;
+    fields          ( U T );
+}
+
+
+// Inflow uniformity, but use a scalar field for weighting
+// Since this field is quite uniform, there should be no difference
+T_UI1
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane1;
+
+    operation       weightedUniformity;
+    weightField     T;
+    fields          ( U );
+}
+
+
+// rhoU-weighted uniformity, including weighting U too (weird but possible)
+rhoU_UI1
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane1;
+
+    operation       weightedUniformity;
+    weightField     rhoU;
+    fields          ( p rho U rhoU );
+}
+
+
+// rhoU-weighted uniformity
+rhoU_UI2
+{
+    ${__surfaceFieldValue}
+
+    regionType      surface;
+    name            plane2;
+
+    operation       weightedUniformity;
+    weightField     rhoU;
+    fields          ( p rho U rhoU );
+}
+
+// ************************************************************************* //