Change of plan - fields will be stored or referenced only at the solver level and tmps will be returned from the clouds to give them their instantaneous values. Only T and U will be supplied to the dsmcCloud, as before, to provide boundary conditions. Added required field calculation functions. Taking care about what is measured - a U field requires sum(momentum)/sum(mass) per cell - cells can have zero particles in them, hence divide by zero problems. Averaging the momentum field and the rhoM field, and constructing the velocity by dividing the averages as a post-processing stage is a safer and more physcially correct method.

2025-11-28 03:28:01 +00:00 · 2009-03-05 19:15:26 +00:00
parent 65bb236e4f
commit a967daeb34
6 changed files with 400 additions and 15 deletions
--- a/applications/solvers/dsmc/dsmcFoam/createFields.H
+++ b/applications/solvers/dsmc/dsmcFoam/createFields.H
@ -1,4 +1,18 @@

+    Info<< nl << "Reading field T" << endl;
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
    Info<< nl << "Reading field U" << endl;
    volVectorField U
    (
@ -13,12 +27,12 @@
        mesh
    );

-    Info<< nl << "Reading field T" << endl;
-    volScalarField T
+    Info<< nl << "Reading field momentum (momentum density)" << endl;
+    volVectorField momentum
    (
        IOobject
        (
-            "T",
+            "momentum",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
@ -104,4 +118,3 @@
    Info<< "Constructing dsmcCloud " << endl;

 dsmcCloud dsmc("dsmc", T, U);
-    //, rhoN, rhoM, dsmcRhoN, q, f);
--- a/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C
+++ b/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C
@ -52,8 +52,26 @@ int main(int argc, char *argv[])

        Info<< "Time = " << runTime.timeName() << nl << endl;

+        // Carry out dsmcCloud timestep
+
        dsmc.evolve();

+        // Retrieve field data from dsmcCloud
+
+        rhoN = dsmc.rhoN();
+
+        rhoM = dsmc.rhoM();
+
+        dsmcRhoN = dsmc.dsmcRhoN();
+
+        momentum = dsmc.momentum();
+
+        q = dsmc.q();
+
+        fD = dsmc.fD();
+
+        // Print status of dsmcCloud
+
        dsmc.info();

        runTime.write();
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
@ -318,6 +318,25 @@ void Foam::DsmcCloud<ParcelType>::collisions()
 }


+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::resetSurfaceDataFields()
+{
+    volScalarField::GeometricBoundaryField& qBF(q_.boundaryField());
+
+    forAll(qBF, p)
+    {
+        qBF[p] = 0.0;
+    }
+
+    volVectorField::GeometricBoundaryField& fDBF(fD_.boundaryField());
+
+    forAll(fDBF, p)
+    {
+        fDBF[p] = vector::zero;
+    }
+}
+
+
 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //

 template<class ParcelType>
@ -385,6 +404,37 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
        mesh_
    ),
    collisionSelectionRemainder_(mesh_.nCells(), 0),
+    q_
+    (
+        IOobject
+        (
+            this->name() + "q_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0)
+    ),
+    fD_
+    (
+        IOobject
+        (
+            this->name() + "fD_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector
+        (
+            "zero",
+            dimensionSet(1, -1, -2, 0, 0),
+            vector::zero
+        )
+    ),
    constProps_(),
    rndGen_(label(971501)),
    T_(T),
@ -459,6 +509,37 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
        dimensionedScalar("zero",  dimensionSet(0, 3, -1, 0, 0), 0.0)
    ),
    collisionSelectionRemainder_(),
+    q_
+    (
+        IOobject
+        (
+            this->name() + "q_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0)
+    ),
+    fD_
+    (
+        IOobject
+        (
+            this->name() + "fD_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector
+        (
+            "zero",
+            dimensionSet(1, -1, -2, 0, 0),
+            vector::zero
+        )
+    ),
    constProps_(),
    rndGen_(label(971501)),
    T_
@ -544,6 +625,9 @@ void Foam::DsmcCloud<ParcelType>::evolve()
    // Insert new particles from the inflow boundary
    this->inflowBoundary().inflow();

+    // Reset the surface data collection fields
+    resetSurfaceDataFields();
+
    // Move the particles ballistically with their current velocities
    Cloud<ParcelType>::move(td);

--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
@ -104,6 +104,12 @@ class DsmcCloud
        //- A field holding the remainder from the previous collision selections
        scalarField collisionSelectionRemainder_;

+        //- Heat flux at surface field
+        volScalarField q_;
+
+        //- Force density at surface field
+        volVectorField fD_;
+
        //- Parcel constant properties - one for each type
        List<typename ParcelType::constantProperties> constProps_;

@ -149,6 +155,9 @@ class DsmcCloud
        //- Calculate collisions between molecules
        void collisions();

+        //- Reset the surface data accumulation field values
+        void resetSurfaceDataFields();
+
        //- Disallow default bitwise copy construct
        DsmcCloud(const DsmcCloud&);

@ -238,6 +247,21 @@ public:
                inline Random& rndGen();


+            // References to the surface data collection fields
+
+                //- Return heat flux at surface field
+                inline const volScalarField& q() const;
+
+                //- Return non-const heat flux at surface field
+                inline volScalarField& q();
+
+                //- Return force density at surface field
+                inline const volVectorField& fD() const;
+
+                //- Return non-const force density at surface field
+                inline volVectorField& fD();
+
+
            // References to the macroscopic fields

                //- Return macroscopic temperature
@ -341,7 +365,20 @@ public:
                inline const tmp<volScalarField> rhoM() const;

                //- Return the field of number of DSMC particles
-                inline const tmp<volScalarField> rhoNP() const;
+                inline const tmp<volScalarField> dsmcRhoN() const;
+
+                //- Return the momentum density field
+                inline const tmp<volVectorField> momentum() const;
+
+                //- Return the total linear kinetic energy (translational and
+                // thermal density field
+                inline const tmp<volScalarField> linearKE() const;
+
+                //- Return the internal energy density field
+                inline const tmp<volScalarField> internalE() const;
+
+                //- Return the average internal degrees of freedom  field
+                inline const tmp<volScalarField> iDof() const;


        // Cloud evolution functions
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
@ -120,6 +120,34 @@ inline Foam::Random& Foam::DsmcCloud<ParcelType>::rndGen()
 }


+template<class ParcelType>
+inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q() const
+{
+    return q_;
+}
+
+
+template<class ParcelType>
+inline Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q()
+{
+    return q_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD() const
+{
+    return fD_;
+}
+
+
+template<class ParcelType>
+inline Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD()
+{
+    return fD_;
+}
+
+
 template<class ParcelType>
 inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::T() const
 {
@ -319,6 +347,17 @@ Foam::DsmcCloud<ParcelType>::rhoN() const
        )
    );

+    scalarField& rhoN = trhoN().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        rhoN[cellI]++;
+    }
+
+    rhoN *= nParticle_/mesh().cellVolumes();
+
    return trhoN;
 }

@ -345,21 +384,32 @@ Foam::DsmcCloud<ParcelType>::rhoM() const
        )
    );

+    scalarField& rhoM = trhoM().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        rhoM[cellI] += constProps(p.typeId()).mass();
+    }
+
+    rhoM *= nParticle_/mesh().cellVolumes();
+
    return trhoM;
 }


 template<class ParcelType>
 inline const Foam::tmp<Foam::volScalarField>
-Foam::DsmcCloud<ParcelType>::rhoNP() const
+Foam::DsmcCloud<ParcelType>::dsmcRhoN() const
 {
-    tmp<volScalarField> trhoNP
+    tmp<volScalarField> tdsmcRhoN
    (
        new volScalarField
        (
            IOobject
            (
-                this->name() + "rhoNP",
+                this->name() + "dsmcRhoN",
                this->db().time().timeName(),
                this->db(),
                IOobject::NO_READ,
@ -371,7 +421,174 @@ Foam::DsmcCloud<ParcelType>::rhoNP() const
        )
    );

-    return trhoNP;
+    scalarField& dsmcRhoN = tdsmcRhoN().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        dsmcRhoN[cellI]++;
+    }
+
+    return tdsmcRhoN;
+}
+
+
+template<class ParcelType>
+inline const Foam::tmp<Foam::volVectorField>
+Foam::DsmcCloud<ParcelType>::momentum() const
+{
+    tmp<volVectorField> tmomentum
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                this->name() + "momentum",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedVector
+            (
+                "zero",
+                dimensionSet(1, -2, -1, 0, 0),
+                vector::zero
+            )
+        )
+    );
+
+    vectorField& momentum = tmomentum().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        momentum[cellI] += constProps(p.typeId()).mass()*p.U();
+    }
+
+    momentum *= nParticle_/mesh().cellVolumes();
+
+    return tmomentum;
+}
+
+
+template<class ParcelType>
+inline const Foam::tmp<Foam::volScalarField>
+Foam::DsmcCloud<ParcelType>::linearKE() const
+{
+    tmp<volScalarField> tlinearKE
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                this->name() + "linearKE",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
+        )
+    );
+
+    scalarField& linearKE = tlinearKE().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        linearKE[cellI] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
+    }
+
+    linearKE *= nParticle_/mesh().cellVolumes();
+
+    return tlinearKE;
+}
+
+
+template<class ParcelType>
+inline const Foam::tmp<Foam::volScalarField>
+Foam::DsmcCloud<ParcelType>::internalE() const
+{
+    tmp<volScalarField> tinternalE
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                this->name() + "internalE",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
+        )
+    );
+
+    scalarField& internalE = tinternalE().internalField();
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        internalE[cellI] += p.Ei();
+    }
+
+    internalE *= nParticle_/mesh().cellVolumes();
+
+    return tinternalE;
+}
+
+
+template<class ParcelType>
+inline const Foam::tmp<Foam::volScalarField>
+Foam::DsmcCloud<ParcelType>::iDof() const
+{
+    tmp<volScalarField> tiDof
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                this->name() + "iDof",
+                this->db().time().timeName(),
+                this->db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimless, 0.0)
+        )
+    );
+
+    scalarField& iDof = tiDof().internalField();
+    scalarField nDsmcParticles(iDof.size(),0);
+
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
+    {
+        const ParcelType& p = iter();
+        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;
+
+    return tiDof;
 }


--- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
@ -109,9 +109,11 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
 {
    const constantProperties& constProps(td.cloud().constProps(typeId_));

-    scalar preInteractionEnergy = 0.5*constProps.mass()*(U_ & U_) + Ei_;
+    // pre-interaction energy
+    scalar preIE = 0.5*constProps.mass()*(U_ & U_) + Ei_;

-    vector preInteractionMomentum = constProps.mass()*U_;
+    // pre-interaction momentum
+    vector preIMom = constProps.mass()*U_;

    td.cloud().wallInteraction().correct
    (
@ -121,13 +123,27 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
        constProps.mass()
    );

-    scalar postInteractionEnergy = 0.5*constProps.mass()*(U_ & U_) + Ei_;
+    // post-interaction energy
+    scalar postIE = 0.5*constProps.mass()*(U_ & U_) + Ei_;

-    vector postInteractionMomentum = constProps.mass()*U_;
+    // post-interaction momentum
+    vector postIMom = constProps.mass()*U_;

-    scalar deltaEnergy =preInteractionEnergy - postInteractionEnergy;
+    label wppIndex = wpp.index();

-    vector deltaMomentum = preInteractionMomentum - postInteractionMomentum;
+    label wppLocalFace = wpp.whichFace(this->face());
+
+    const scalar fA = mag(wpp.faceAreas()[wppLocalFace]);
+
+    const scalar deltaT = td.cloud().mesh().time().deltaT().value();
+
+    scalar deltaQ = td.cloud().nParticle()*(preIE - postIE)/(deltaT*fA);
+
+    vector deltaFD = td.cloud().nParticle()*(preIMom - postIMom)/(deltaT*fA);
+
+    td.cloud().q().boundaryField()[wppIndex][wppLocalFace] += deltaQ;
+
+    td.cloud().fD().boundaryField()[wppIndex][wppLocalFace] += deltaFD;
 }