diff --git a/applications/solvers/dsmc/dsmcFoam/createFields.H b/applications/solvers/dsmc/dsmcFoam/createFields.H index 5e5751bddf..74c1d5c055 100644 --- 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); diff --git a/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C b/applications/solvers/dsmc/dsmcFoam/dsmcFoam.C index 91b53b43a9..a9d30e9ff0 100644 --- 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(); diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C index c588bb0cb7..0e39c2ff33 100644 --- 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::collisions() } +template +void Foam::DsmcCloud::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 @@ -385,6 +404,37 @@ Foam::DsmcCloud::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::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::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::move(td); diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H index d0bb26bc9c..8e8fd5cbd7 100644 --- 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 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 rhoM() const; //- Return the field of number of DSMC particles - inline const tmp rhoNP() const; + inline const tmp dsmcRhoN() const; + + //- Return the momentum density field + inline const tmp momentum() const; + + //- Return the total linear kinetic energy (translational and + // thermal density field + inline const tmp linearKE() const; + + //- Return the internal energy density field + inline const tmp internalE() const; + + //- Return the average internal degrees of freedom field + inline const tmp iDof() const; // Cloud evolution functions diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H index 39c25f5227..d5dda1f680 100644 --- 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::rndGen() } +template +inline const Foam::volScalarField& Foam::DsmcCloud::q() const +{ + return q_; +} + + +template +inline Foam::volScalarField& Foam::DsmcCloud::q() +{ + return q_; +} + + +template +inline const Foam::volVectorField& Foam::DsmcCloud::fD() const +{ + return fD_; +} + + +template +inline Foam::volVectorField& Foam::DsmcCloud::fD() +{ + return fD_; +} + + template inline const Foam::volScalarField& Foam::DsmcCloud::T() const { @@ -319,6 +347,17 @@ Foam::DsmcCloud::rhoN() const ) ); + scalarField& rhoN = trhoN().internalField(); + forAllConstIter(typename DsmcCloud, *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::rhoM() const ) ); + scalarField& rhoM = trhoM().internalField(); + forAllConstIter(typename DsmcCloud, *this, iter) + { + const ParcelType& p = iter(); + const label cellI = p.cell(); + + rhoM[cellI] += constProps(p.typeId()).mass(); + } + + rhoM *= nParticle_/mesh().cellVolumes(); + return trhoM; } template inline const Foam::tmp -Foam::DsmcCloud::rhoNP() const +Foam::DsmcCloud::dsmcRhoN() const { - tmp trhoNP + tmp 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::rhoNP() const ) ); - return trhoNP; + scalarField& dsmcRhoN = tdsmcRhoN().internalField(); + forAllConstIter(typename DsmcCloud, *this, iter) + { + const ParcelType& p = iter(); + const label cellI = p.cell(); + + dsmcRhoN[cellI]++; + } + + return tdsmcRhoN; +} + + +template +inline const Foam::tmp +Foam::DsmcCloud::momentum() const +{ + tmp 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, *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 +inline const Foam::tmp +Foam::DsmcCloud::linearKE() const +{ + tmp 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, *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 +inline const Foam::tmp +Foam::DsmcCloud::internalE() const +{ + tmp 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, *this, iter) + { + const ParcelType& p = iter(); + const label cellI = p.cell(); + + internalE[cellI] += p.Ei(); + } + + internalE *= nParticle_/mesh().cellVolumes(); + + return tinternalE; +} + + +template +inline const Foam::tmp +Foam::DsmcCloud::iDof() const +{ + tmp 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, *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; } diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C index babf4be292..93058e076e 100644 --- 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::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::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; }