Averaging does not need to be weighted if the fields being averaged are extensive - they weight themselves automatically. Recording momentum and energy denisty, and averaging them, then once decent averaged fields exist, using the new dsmcFields utility to create the intensive fields: UMean and TMean (translational, internal and overall).

applications/solvers/dsmc/dsmcFoam/createFields.H

@@ -1,10 +1,10 @@
 
-    Info<< nl << "Reading field T" << endl;
+    Info<< nl << "Reading field boundaryT" << endl;
-    volScalarField T
+    volScalarField boundaryT
     (
         IOobject
         (
-            "T",
+            "boundaryT",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -13,12 +13,12 @@
         mesh
     );
 
-    Info<< nl << "Reading field U" << endl;
+    Info<< nl << "Reading field boundaryU" << endl;
-    volVectorField U
+    volVectorField boundaryU
     (
         IOobject
         (
-            "U",
+            "boundaryU",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -27,21 +27,7 @@
         mesh
     );
 
-    Info<< nl << "Reading field momentum (momentum density)" << endl;
+    Info<< nl << "Reading field rhoN (number density)" << endl;
-    volVectorField momentum
-    (
-        IOobject
-        (
-            "momentum",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Reading field rhoN" << endl;
     volScalarField rhoN
     (
         IOobject
@@ -55,7 +41,7 @@
         mesh
     );
 
-    Info<< nl << "Reading field rhoM" << endl;
+    Info<< nl << "Reading field rhoM (mass density)" << endl;
     volScalarField rhoM
     (
         IOobject
@@ -83,10 +69,66 @@
         mesh
     );
 
-    Info<< nl
+    Info<< nl << "Reading field momentum (momentum density)" << endl;
-        << "Fields required to record solid surface forces and heat flux:"
+    volVectorField momentum
+    (
+        IOobject
+        (
+            "momentum",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field linearKE (linear kinetic energy density)"
         << endl;
 
+    volScalarField linearKE
+    (
+        IOobject
+        (
+            "linearKE",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field internalE (internal energy density)" << endl;
+    volScalarField internalE
+    (
+        IOobject
+        (
+            "internalE",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< nl << "Reading field iDof (internal degree of freedom density)"
+        << endl;
+
+    volScalarField iDof
+    (
+        IOobject
+        (
+            "iDof",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
     Info<< nl << "Reading field q (surface heat transfer)" << endl;
     volScalarField q
     (
@@ -101,7 +143,7 @@
         mesh
     );
 
-    Info<< nl << "Reading field fD (surface force density" << endl;
+    Info<< nl << "Reading field fD (surface force density)" << endl;
     volVectorField fD
     (
         IOobject
@@ -115,6 +157,6 @@
         mesh
     );
 
-    Info<< "Constructing dsmcCloud " << endl;
+    Info<< nl << "Constructing dsmcCloud " << endl;
 
-dsmcCloud dsmc("dsmc", T, U);
+    dsmcCloud dsmc("dsmc", boundaryT, boundaryU);

applications/solvers/dsmc/dsmcFoam/dsmcFoam.C

@@ -66,6 +66,12 @@ int main(int argc, char *argv[])
 
         momentum = dsmc.momentum();
 
+        linearKE = dsmc.linearKE();
+
+        internalE = dsmc.internalE();
+
+        iDof = dsmc.iDof();
+
         q = dsmc.q();
 
         fD = dsmc.fD();

applications/utilities/postProcessing/miscellaneous/dsmcFields/Make/files

@@ -0,0 +1,3 @@
+dsmcFields.C
+
+EXE = $(FOAM_APPBIN)/dsmcFields

applications/utilities/postProcessing/miscellaneous/dsmcFields/Make/options

@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(LIB_SRC)/postProcessing/postCalc \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/dsmc/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    $(FOAM_LIBBIN)/postCalc.o \
+    -lmeshTools \
+    -lfiniteVolume \
+    -llagrangian \
+    -ldsmc
+

applications/utilities/postProcessing/miscellaneous/dsmcFields/dsmcFields.C

@@ -0,0 +1,231 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    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
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    dsmcFields
+
+Description
+    Calculate intensive fields (U and T) from averaged extensive fields from a
+    DSMC calculation.
+
+\*---------------------------------------------------------------------------*/
+
+#include "calc.H"
+#include "fvc.H"
+#include "dsmcCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
+{
+    bool writeResults = !args.options().found("noWrite");
+
+    IOobject rhoNMeanheader
+    (
+        "rhoNMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    IOobject rhoMMeanheader
+    (
+        "rhoMMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    IOobject momentumMeanheader
+    (
+        "momentumMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    IOobject linearKEMeanheader
+    (
+        "linearKEMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    IOobject internalEMeanheader
+    (
+        "internalEMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    IOobject iDofMeanheader
+    (
+        "iDofMean",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ
+    );
+
+    if (!rhoNMeanheader.headerOk())
+    {
+        Info<< "    No rhoNMean" << endl;
+    }
+    else if (!rhoMMeanheader.headerOk())
+    {
+        Info<< "    No rhoMMean" << endl;
+    }
+    else if (!momentumMeanheader.headerOk())
+    {
+        Info<< "    No momentumMean" << endl;
+    }
+    else if (!linearKEMeanheader.headerOk())
+    {
+        Info<< "    No linearKEMean" << endl;
+    }
+    else if (!internalEMeanheader.headerOk())
+    {
+        Info<< "    No internalEMean" << endl;
+    }
+    else if (!iDofMeanheader.headerOk())
+    {
+        Info<< "    No iDofMean" << endl;
+    }
+    else
+    {
+        Info<< "Reading field rhoNMean" << endl;
+        volScalarField rhoNMean(rhoNMeanheader, mesh);
+
+        Info<< "Reading field rhoMMean" << endl;
+        volScalarField rhoMMean(rhoMMeanheader, mesh);
+
+        Info<< "Reading field momentumMean" << endl;
+        volVectorField momentumMean(momentumMeanheader, mesh);
+
+        Info<< "Reading field linearKEMean" << endl;
+        volScalarField linearKEMean(linearKEMeanheader, mesh);
+
+        Info<< "Reading field internalEMean" << endl;
+        volScalarField internalEMean(internalEMeanheader, mesh);
+
+        Info<< "Reading field iDofMean" << endl;
+        volScalarField iDofMean(iDofMeanheader, mesh);
+
+        // Check if there are any zero values in the density fields
+
+        if (min(rhoNMean).value() > VSMALL)
+        {
+            // TODO Sort out boundary field values if required
+
+            Info<< nl << "Calculating UMean field." << endl;
+
+            volVectorField UMean
+            (
+                IOobject
+                (
+                    "UMean",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ
+                ),
+                momentumMean/rhoMMean
+            );
+
+            Info<< nl << "Calculating TMean fields." << endl;
+
+            volScalarField translationalTMean
+            (
+                IOobject
+                (
+                    "translationalTMean",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ
+                ),
+                2.0/(3.0*dsmcCloud::kb*rhoNMean)
+               *(linearKEMean - 0.5*rhoMMean*(UMean & UMean))
+            );
+
+            volScalarField internalTMean
+            (
+                IOobject
+                (
+                    "internalTMean",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ
+                ),
+                2.0/(dsmcCloud::kb*iDofMean)*internalEMean
+            );
+
+            volScalarField overallTMean
+            (
+                IOobject
+                (
+                    "overallTMean",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ
+                ),
+                2.0/(dsmcCloud::kb*(3.0*rhoNMean + iDofMean))
+               *(linearKEMean - 0.5*rhoMMean*(UMean & UMean) + internalEMean)
+            );
+
+            Info<< nl << "magUMean max/min : "
+                << max(mag(UMean)).value() << " "
+                << min(mag(UMean)).value() << endl;
+
+            Info<< nl << "translationalTMean max/min : "
+                << max(translationalTMean).value() << " "
+                << min(translationalTMean).value() << endl;
+
+            Info<< nl << "internalTMean max/min : "
+                << max(internalTMean).value() << " "
+                << min(internalTMean).value() << endl;
+
+            Info<< nl << "overallTMean max/min : "
+                << max(overallTMean).value() << " "
+                << min(overallTMean).value() << endl;
+
+            if (writeResults)
+            {
+                UMean.write();
+                translationalTMean.write();
+                internalTMean.write();
+                overallTMean.write();
+            }
+        }
+        else
+        {
+            Info<< "Small or negative value (" << min(rhoNMean)
+                << ") found in rhoNMean field. "
+                << "Not calculating fields to avoid division by zero "
+                << "or invalid results."
+                << endl;
+        }
+    }
+}
+
+// ************************************************************************* //

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H

@@ -568,12 +568,11 @@ Foam::DsmcCloud<ParcelType>::iDof() const
                 false
             ),
             mesh_,
-            dimensionedScalar("zero", dimless, 0.0)
+            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
         )
     );
 
     scalarField& iDof = tiDof().internalField();
-    scalarField nDsmcParticles(iDof.size(),0);
 
     forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
     {
@@ -581,12 +580,9 @@ Foam::DsmcCloud<ParcelType>::iDof() const
         const label cellI = p.cell();
 
         iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
-
-        // Avoiding divide by zero for nDsmcParticles where a cell is empty
-        nDsmcParticles[cellI] = max(1, ++nDsmcParticles[cellI]);
     }
 
-    iDof /= nDsmcParticles;
+    iDof *= nParticle_/mesh().cellVolumes();
 
     return tiDof;
 }