From cadba3c278e9e581d9aa4c10ef6cc7290c046819 Mon Sep 17 00:00:00 2001 From: graham Date: Tue, 14 Jul 2009 18:46:46 +0100 Subject: [PATCH 01/16] Adding new MixedDiffuseSpecular wall collision model. --- .../makeDsmcParcelWallInteractionModels.C | 7 + .../MixedDiffuseSpecular.C | 140 ++++++++++++++++++ .../MixedDiffuseSpecular.H | 106 +++++++++++++ 3 files changed, 253 insertions(+) create mode 100644 src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C create mode 100644 src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C index f15d0bdadd..b67f66ec85 100644 --- a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C +++ b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C @@ -28,6 +28,7 @@ License #include "DsmcCloud.H" #include "MaxwellianThermal.H" #include "SpecularReflection.H" +#include "MixedDiffuseSpecular.H" namespace Foam { @@ -46,6 +47,12 @@ namespace Foam DsmcCloud, dsmcParcel ); + makeWallInteractionModelType + ( + MixedDiffuseSpecular, + DsmcCloud, + dsmcParcel + ); }; diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C new file mode 100644 index 0000000000..a12fff43d1 --- /dev/null +++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C @@ -0,0 +1,140 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2009-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 + +\*---------------------------------------------------------------------------*/ + +#include "MixedDiffuseSpecular.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::MixedDiffuseSpecular::MixedDiffuseSpecular +( + const dictionary& dict, + CloudType& cloud +) +: + WallInteractionModel(dict, cloud, typeName), + diffuseFraction_(readScalar(this->coeffDict().lookup("diffuseFraction"))) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::MixedDiffuseSpecular::~MixedDiffuseSpecular() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::MixedDiffuseSpecular::correct +( + const wallPolyPatch& wpp, + const label faceId, + vector& U, + scalar& Ei, + label typeId +) +{ + label wppIndex = wpp.index(); + + label wppLocalFace = wpp.whichFace(faceId); + + vector nw = wpp.faceAreas()[wppLocalFace]; + + // Normal unit vector + nw /= mag(nw); + + // Normal velocity magnitude + scalar magUn = U & nw; + + CloudType& cloud(this->owner()); + + Random& rndGen(cloud.rndGen()); + + if (diffuseFraction_ > rndGen.scalar01()) + { + // Diffuse reflection + + // Wall tangential velocity (flow direction) + vector Ut = U - magUn*nw; + + while (mag(Ut) < SMALL) + { + // If the incident velocity is parallel to the face normal, no + // tangential direction can be chosen. Add a perturbation to the + // incoming velocity and recalculate. + + U = vector + ( + U.x()*(0.8 + 0.2*rndGen.scalar01()), + U.y()*(0.8 + 0.2*rndGen.scalar01()), + U.z()*(0.8 + 0.2*rndGen.scalar01()) + ); + + magUn = U & nw; + + Ut = U - magUn*nw; + } + + // Wall tangential unit vector + vector tw1 = Ut/mag(Ut); + + // Other tangential unit vector + vector tw2 = nw^tw1; + + scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace]; + + scalar mass = cloud.constProps(typeId).mass(); + + scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom(); + + U = + sqrt(CloudType::kb*T/mass) + *( + rndGen.GaussNormal()*tw1 + + rndGen.GaussNormal()*tw2 + - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw + ); + + U += cloud.U().boundaryField()[wppIndex][wppLocalFace]; + + Ei = cloud.equipartitionInternalEnergy(T, iDof); + } + else + { + // Specular reflection + + if (magUn > 0.0) + { + U -= 2.0*magUn*nw; + } + } + +} + + +// ************************************************************************* // diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H new file mode 100644 index 0000000000..b84f43d8f6 --- /dev/null +++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H @@ -0,0 +1,106 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2009-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 + +Class + Foam::MixedDiffuseSpecular + +Description + Wall interaction setting microscopic velocity to a random one drawn from a + Maxwellian distribution corresponding to a specified temperature + +\*---------------------------------------------------------------------------*/ + +#ifndef MixedDiffuseSpecular_H +#define MixedDiffuseSpecular_H + +#include "WallInteractionModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +/*---------------------------------------------------------------------------*\ + Class MixedDiffuseSpecular Declaration +\*---------------------------------------------------------------------------*/ + +template +class MixedDiffuseSpecular +: + public WallInteractionModel +{ + // Private data + + //- Fraction of wall interactions that are diffuse + scalar diffuseFraction_; + + +public: + + //- Runtime type information + TypeName("MixedDiffuseSpecular"); + + + // Constructors + + //- Construct from dictionary + MixedDiffuseSpecular + ( + const dictionary& dict, + CloudType& cloud + ); + + + // Destructor + virtual ~MixedDiffuseSpecular(); + + + // Member Functions + + //- Apply wall correction + virtual void correct + ( + const wallPolyPatch& wpp, + const label faceId, + vector& U, + scalar& Ei, + label typeId + ); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "MixedDiffuseSpecular.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // From 9321f7e1e5bf1ee919cf16a7f1e7f388d6f7133e Mon Sep 17 00:00:00 2001 From: graham Date: Wed, 15 Jul 2009 15:28:04 +0100 Subject: [PATCH 02/16] Adding pressure field measurement, internal and surface. --- .../dsmcFieldsCalc/dsmcFieldsCalc.C | 11 ++++- .../utilities/dsmcFields/dsmcFields.C | 43 ++++++++++++++++++- 2 files changed, 51 insertions(+), 3 deletions(-) diff --git a/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C b/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C index 95077df225..34f231850f 100644 --- a/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C +++ b/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C @@ -113,7 +113,16 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh) return; } - wordList extensiveVVFNames(IStringStream ("(momentumMean)")()); + wordList extensiveVVFNames + ( + IStringStream + ( + "( \ + momentumMean \ + fDMean \ + )" + )() + ); PtrList extensiveVVFs(extensiveVVFNames.size()); diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C index 8957dee17c..7056fba521 100644 --- a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C +++ b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C @@ -106,6 +106,7 @@ void Foam::dsmcFields::write() word linearKEMeanName = "linearKEMean"; word internalEMeanName = "internalEMean"; word iDofMeanName = "iDofMean"; + word fDMeanName = "fDMean"; const volScalarField& rhoNMean = obr_.lookupObject ( @@ -137,6 +138,11 @@ void Foam::dsmcFields::write() iDofMeanName ); + volVectorField fDMean = obr_.lookupObject + ( + fDMeanName + ); + if (min(mag(rhoNMean)).value() > VSMALL) { Info<< "Calculating dsmcFields." << endl; @@ -165,7 +171,7 @@ void Foam::dsmcFields::write() IOobject::NO_READ ), 2.0/(3.0*dsmcCloud::kb*rhoNMean) - *(linearKEMean - 0.5*rhoMMean*(UMean & UMean)) + *(linearKEMean - 0.5*rhoMMean*(UMean & UMean)) ); Info<< " Calculating internalT field." << endl; @@ -192,9 +198,36 @@ void Foam::dsmcFields::write() IOobject::NO_READ ), 2.0/(dsmcCloud::kb*(3.0*rhoNMean + iDofMean)) - *(linearKEMean - 0.5*rhoMMean*(UMean & UMean) + internalEMean) + *(linearKEMean - 0.5*rhoMMean*(UMean & UMean) + internalEMean) ); + Info<< " Calculating pressure field." << endl; + volScalarField p + ( + IOobject + ( + "p", + obr_.time().timeName(), + obr_, + IOobject::NO_READ + ), + dsmcCloud::kb*rhoNMean*translationalT + ); + + const fvMesh& mesh = fDMean.mesh(); + + forAll(mesh.boundaryMesh(), i) + { + const polyPatch& patch = mesh.boundaryMesh()[i]; + + if (isA(patch)) + { + p.boundaryField()[i] = + fDMean.boundaryField()[i] + & (patch.faceAreas()/mag(patch.faceAreas())); + } + } + Info<< " mag(UMean) max/min : " << max(mag(UMean)).value() << " " << min(mag(UMean)).value() << endl; @@ -211,6 +244,10 @@ void Foam::dsmcFields::write() << max(overallT).value() << " " << min(overallT).value() << endl; + Info<< " p max/min : " + << max(p).value() << " " + << min(p).value() << endl; + UMean.write(); translationalT.write(); @@ -219,6 +256,8 @@ void Foam::dsmcFields::write() overallT.write(); + p.write(); + Info<< "dsmcFields written." << nl << endl; } else From 103cd27ef563091d959268b9cc92691662e3dfb7 Mon Sep 17 00:00:00 2001 From: graham Date: Fri, 17 Jul 2009 15:40:11 +0100 Subject: [PATCH 03/16] Measurement of surface velocity. Created momentum and rhoM field inside the cloud - duplicated memory with solver for the moment. --- .../clouds/Templates/DsmcCloud/DsmcCloud.C | 79 ++++++++ .../clouds/Templates/DsmcCloud/DsmcCloud.H | 40 +++-- .../clouds/Templates/DsmcCloud/DsmcCloudI.H | 169 ++++++++---------- .../parcels/Templates/DsmcParcel/DsmcParcel.C | 42 +++-- .../MaxwellianThermal/MaxwellianThermal.C | 8 +- .../MixedDiffuseSpecular.C | 12 +- .../SpecularReflection/SpecularReflection.C | 6 +- 7 files changed, 225 insertions(+), 131 deletions(-) diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C index 23aa3dc478..8615135a26 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C @@ -490,6 +490,23 @@ void Foam::DsmcCloud::resetSurfaceDataFields() { fDBF[p] = vector::zero; } + + volVectorField::GeometricBoundaryField& momentumBF + ( + momentum_.boundaryField() + ); + + forAll(momentumBF, p) + { + momentumBF[p] = vector::zero; + } + + volScalarField::GeometricBoundaryField& rhoMBF(rhoM_.boundaryField()); + + forAll(rhoMBF, p) + { + rhoMBF[p] = VSMALL; + } } @@ -591,6 +608,37 @@ Foam::DsmcCloud::DsmcCloud vector::zero ) ), + momentum_ + ( + IOobject + ( + this->name() + "momentum_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedVector + ( + "zero", + dimensionSet(1, -2, -1, 0, 0), + vector::zero + ) + ), + rhoM_ + ( + IOobject + ( + this->name() + "rhoM_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), 0.0) + ), constProps_(), rndGen_(label(149382906) + 7183*Pstream::myProcNo()), T_(T), @@ -703,6 +751,37 @@ Foam::DsmcCloud::DsmcCloud vector::zero ) ), + momentum_ + ( + IOobject + ( + this->name() + "momentum_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedVector + ( + "zero", + dimensionSet(1, -2, -1, 0, 0), + vector::zero + ) + ), + rhoM_ + ( + IOobject + ( + this->name() + "rhoM_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), 0.0) + ), constProps_(), rndGen_(label(971501) + 1526*Pstream::myProcNo()), T_ diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H index 76b1dd7e79..40f3935eee 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H @@ -110,6 +110,12 @@ class DsmcCloud //- Force density at surface field volVectorField fD_; + //- Momentum density field + volVectorField momentum_; + + //- Mass density field + volScalarField rhoM_; + //- Parcel constant properties - one for each type List constProps_; @@ -254,19 +260,19 @@ public: inline scalar cachedDeltaT() const; - // References to the surface data collection fields + // References to the boundary fields for surface data collection - //- Return heat flux at surface field - inline const volScalarField& q() const; + //- Return non-const heat flux boundary field reference + inline volScalarField::GeometricBoundaryField& qBF(); - //- Return non-const heat flux at surface field - inline volScalarField& q(); + //- Return non-const force density at boundary field reference + inline volVectorField::GeometricBoundaryField& fDBF(); - //- Return force density at surface field - inline const volVectorField& fD() const; + //- Return non-const momentum density boundary field reference + inline volVectorField::GeometricBoundaryField& momentumBF(); - //- Return non-const force density at surface field - inline volVectorField& fD(); + //- Return non-const mass density boundary field reference + inline volScalarField::GeometricBoundaryField& rhoMBF(); // References to the macroscopic fields @@ -391,17 +397,20 @@ public: // Fields + //- Return heat flux at surface field + inline const volScalarField& q() const; + + //- Return force density at surface field + inline const volVectorField& fD() const; + //- Return the real particle number density field inline const tmp rhoN() const; //- Return the particle mass density field - inline const tmp rhoM() const; - - //- Return the field of number of DSMC particles - inline const tmp dsmcRhoN() const; + inline const volScalarField& rhoM(); //- Return the momentum density field - inline const tmp momentum() const; + inline const volVectorField& momentum(); //- Return the total linear kinetic energy (translational and // thermal density field @@ -413,6 +422,9 @@ public: //- Return the average internal degrees of freedom field inline const tmp iDof() const; + //- Return the field of number of DSMC particles + inline const tmp dsmcRhoN() 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 3529702fdd..07ac87bd2c 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H @@ -135,30 +135,34 @@ inline Foam::scalar Foam::DsmcCloud::cachedDeltaT() const template -inline const Foam::volScalarField& Foam::DsmcCloud::q() const +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::qBF() { - return q_; + return q_.boundaryField(); } template -inline Foam::volScalarField& Foam::DsmcCloud::q() +inline Foam::volVectorField::GeometricBoundaryField& +Foam::DsmcCloud::fDBF() { - return q_; + return fD_.boundaryField(); } template -inline const Foam::volVectorField& Foam::DsmcCloud::fD() const +inline Foam::volVectorField::GeometricBoundaryField& +Foam::DsmcCloud::momentumBF() { - return fD_; + return momentum_.boundaryField(); } template -inline Foam::volVectorField& Foam::DsmcCloud::fD() +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::rhoMBF() { - return fD_; + return rhoM_.boundaryField(); } @@ -373,6 +377,20 @@ Foam::DsmcCloud::maxwellianMostProbableSpeed } +template +inline const Foam::volScalarField& Foam::DsmcCloud::q() const +{ + return q_; +} + + +template +inline const Foam::volVectorField& Foam::DsmcCloud::fD() const +{ + return fD_; +} + + template inline const Foam::tmp Foam::DsmcCloud::rhoN() const @@ -411,116 +429,44 @@ Foam::DsmcCloud::rhoN() const template -inline const Foam::tmp -Foam::DsmcCloud::rhoM() const +inline const Foam::volScalarField& Foam::DsmcCloud::rhoM() { - tmp trhoM - ( - new volScalarField - ( - IOobject - ( - this->name() + "rhoM", - this->db().time().timeName(), - this->db(), - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh_, - dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), VSMALL) - ) - ); + scalarField& rM = rhoM_.internalField(); + + rM = VSMALL; - scalarField& rhoM = trhoM().internalField(); forAllConstIter(typename DsmcCloud, *this, iter) { const ParcelType& p = iter(); const label cellI = p.cell(); - rhoM[cellI] += constProps(p.typeId()).mass(); + rM[cellI] += constProps(p.typeId()).mass(); } - rhoM *= nParticle_/mesh().cellVolumes(); + rM *= nParticle_/mesh().cellVolumes(); - return trhoM; + return rhoM_; } template -inline const Foam::tmp -Foam::DsmcCloud::dsmcRhoN() const +inline const Foam::volVectorField& Foam::DsmcCloud::momentum() { - tmp tdsmcRhoN - ( - new volScalarField - ( - IOobject - ( - this->name() + "dsmcRhoN", - this->db().time().timeName(), - this->db(), - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh_, - dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0) - ) - ); + vectorField& mom = momentum_.internalField(); + + mom = vector::zero; - scalarField& dsmcRhoN = tdsmcRhoN().internalField(); forAllConstIter(typename DsmcCloud, *this, iter) { const ParcelType& p = iter(); const label cellI = p.cell(); - dsmcRhoN[cellI]++; + mom[cellI] += constProps(p.typeId()).mass()*p.U(); } - return tdsmcRhoN; -} + mom *= nParticle_/mesh().cellVolumes(); - -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; + return momentum_; } @@ -636,6 +582,41 @@ Foam::DsmcCloud::iDof() const } +template +inline const Foam::tmp +Foam::DsmcCloud::dsmcRhoN() const +{ + tmp tdsmcRhoN + ( + new volScalarField + ( + IOobject + ( + this->name() + "dsmcRhoN", + this->db().time().timeName(), + this->db(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh_, + dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0) + ) + ); + + scalarField& dsmcRhoN = tdsmcRhoN().internalField(); + forAllConstIter(typename DsmcCloud, *this, iter) + { + const ParcelType& p = iter(); + const label cellI = p.cell(); + + dsmcRhoN[cellI]++; + } + + return tdsmcRhoN; +} + + template inline void Foam::DsmcCloud::clear() { diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C index c52d92d20f..e11a5e73ea 100644 --- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C +++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C @@ -114,10 +114,30 @@ void Foam::DsmcParcel::hitWallPatch TrackData& td ) { + label wppIndex = wpp.index(); + + label wppLocalFace = wpp.whichFace(this->face()); + + const scalar fA = mag(wpp.faceAreas()[wppLocalFace]); + + const scalar deltaT = td.cloud().cachedDeltaT(); + const constantProperties& constProps(td.cloud().constProps(typeId_)); scalar m = constProps.mass(); + vector nw = wpp.faceAreas()[wppLocalFace]; + nw /= mag(nw); + + scalar U_dot_nw = U_ & nw; + vector Ut = U_ - U_dot_nw*nw; + + td.cloud().momentumBF()[wppIndex][wppLocalFace] += + m*Ut/max(mag(U_dot_nw)*fA, VSMALL); + + td.cloud().rhoMBF()[wppIndex][wppLocalFace] += + m/max(mag(U_dot_nw)*fA, VSMALL); + // pre-interaction energy scalar preIE = 0.5*m*(U_ & U_) + Ei_; @@ -133,27 +153,29 @@ void Foam::DsmcParcel::hitWallPatch typeId_ ); + U_dot_nw = U_ & nw; + Ut = U_ - U_dot_nw*nw; + + td.cloud().momentumBF()[wppIndex][wppLocalFace] += + m*Ut/max(mag(U_dot_nw)*fA, VSMALL); + + td.cloud().rhoMBF()[wppIndex][wppLocalFace] += + m/max(mag(U_dot_nw)*fA, VSMALL); + // post-interaction energy scalar postIE = 0.5*m*(U_ & U_) + Ei_; // post-interaction momentum vector postIMom = m*U_; - label wppIndex = wpp.index(); - - label wppLocalFace = wpp.whichFace(this->face()); - - const scalar fA = mag(wpp.faceAreas()[wppLocalFace]); - - const scalar deltaT = td.cloud().cachedDeltaT(); - 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().qBF()[wppIndex][wppLocalFace] += deltaQ; + + td.cloud().fDBF()[wppIndex][wppLocalFace] += deltaFD; - td.cloud().fD().boundaryField()[wppIndex][wppLocalFace] += deltaFD; } diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C index e7213561c0..6ad452793b 100644 --- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C +++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C @@ -68,10 +68,10 @@ void Foam::MaxwellianThermal::correct nw /= mag(nw); // Normal velocity magnitude - scalar magUn = U & nw; + scalar U_dot_nw = U & nw; // Wall tangential velocity (flow direction) - vector Ut = U - magUn*nw; + vector Ut = U - U_dot_nw*nw; CloudType& cloud(this->owner()); @@ -90,9 +90,9 @@ void Foam::MaxwellianThermal::correct U.z()*(0.8 + 0.2*rndGen.scalar01()) ); - magUn = U & nw; + U_dot_nw = U & nw; - Ut = U - magUn*nw; + Ut = U - U_dot_nw*nw; } // Wall tangential unit vector diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C index a12fff43d1..d491dfa406 100644 --- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C +++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C @@ -69,7 +69,7 @@ void Foam::MixedDiffuseSpecular::correct nw /= mag(nw); // Normal velocity magnitude - scalar magUn = U & nw; + scalar U_dot_nw = U & nw; CloudType& cloud(this->owner()); @@ -80,7 +80,7 @@ void Foam::MixedDiffuseSpecular::correct // Diffuse reflection // Wall tangential velocity (flow direction) - vector Ut = U - magUn*nw; + vector Ut = U - U_dot_nw*nw; while (mag(Ut) < SMALL) { @@ -95,9 +95,9 @@ void Foam::MixedDiffuseSpecular::correct U.z()*(0.8 + 0.2*rndGen.scalar01()) ); - magUn = U & nw; + U_dot_nw = U & nw; - Ut = U - magUn*nw; + Ut = U - U_dot_nw*nw; } // Wall tangential unit vector @@ -128,9 +128,9 @@ void Foam::MixedDiffuseSpecular::correct { // Specular reflection - if (magUn > 0.0) + if (U_dot_nw > 0.0) { - U -= 2.0*magUn*nw; + U -= 2.0*U_dot_nw*nw; } } diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C index 1fb5ab3208..bed84b8c68 100644 --- a/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C +++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C @@ -63,11 +63,11 @@ void Foam::SpecularReflection::correct vector nw = wpp.faceAreas()[wpp.whichFace(faceId)]; nw /= mag(nw); - scalar magUn = U & nw; + scalar U_dot_nw = U & nw; - if (magUn > 0.0) + if (U_dot_nw > 0.0) { - U -= 2.0*magUn*nw; + U -= 2.0*U_dot_nw*nw; } } From 85d7852fc350270c175090e8ccd1e45b4b99b056 Mon Sep 17 00:00:00 2001 From: graham Date: Fri, 17 Jul 2009 18:17:26 +0100 Subject: [PATCH 04/16] Renamed U_ and T_ to boundaryU_ and boundaryT_. Moved all fields and field reading into the DsmcCloud, all calculation and resetting to single functions for all fields. Changed constructors so that no fields are supplied to the solver called from dsmcFoam and an initialisation dictionary is supplied by dsmcInitialise. --- .../dsmc/dsmcFoam/createFields.H | 162 -------- .../discreteMethods/dsmc/dsmcFoam/dsmcFoam.C | 42 +- .../dsmcInitialise/dsmcInitialise.C | 14 +- .../clouds/Templates/DsmcCloud/DsmcCloud.C | 375 +++++++++++++----- .../clouds/Templates/DsmcCloud/DsmcCloud.H | 123 +++--- .../clouds/Templates/DsmcCloud/DsmcCloudI.H | 248 ++---------- .../dsmc/clouds/derived/dsmcCloud/dsmcCloud.C | 10 +- .../dsmc/clouds/derived/dsmcCloud/dsmcCloud.H | 19 +- .../FreeStream/FreeStream.C | 7 +- .../MaxwellianThermal/MaxwellianThermal.C | 4 +- .../MixedDiffuseSpecular.C | 4 +- 11 files changed, 433 insertions(+), 575 deletions(-) delete mode 100644 applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H diff --git a/applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H b/applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H deleted file mode 100644 index d024bd2017..0000000000 --- a/applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H +++ /dev/null @@ -1,162 +0,0 @@ - - Info<< nl << "Reading field boundaryT" << endl; - volScalarField boundaryT - ( - IOobject - ( - "boundaryT", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - Info<< nl << "Reading field boundaryU" << endl; - volVectorField boundaryU - ( - IOobject - ( - "boundaryU", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - Info<< nl << "Reading field rhoN (number density)" << endl; - volScalarField rhoN - ( - IOobject - ( - "rhoN", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - Info<< nl << "Reading field rhoM (mass density)" << endl; - volScalarField rhoM - ( - IOobject - ( - "rhoM", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - Info<< nl << "Reading field rhoNdsmc (dsmc particle density)" << endl; - volScalarField dsmcRhoN - ( - IOobject - ( - "dsmcRhoN", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - Info<< nl << "Reading field momentum (momentum density)" << endl; - 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 - ( - IOobject - ( - "q", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - Info<< nl << "Reading field fD (surface force density)" << endl; - volVectorField fD - ( - IOobject - ( - "fD", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - Info<< nl << "Constructing dsmcCloud " << endl; - - dsmcCloud dsmc("dsmc", boundaryT, boundaryU); diff --git a/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C b/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C index 7a4d8ce8b9..89c7c817e9 100644 --- a/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C +++ b/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C @@ -40,53 +40,21 @@ int main(int argc, char *argv[]) #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" - #include "createFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + Info<< nl << "Constructing dsmcCloud " << endl; + + dsmcCloud dsmc("dsmc", mesh); + Info<< "\nStarting time loop\n" << endl; - while (runTime.run()) + while (runTime.loop()) { - runTime++; - Info<< "Time = " << runTime.timeName() << nl << endl; - // Carry out dsmcCloud timestep - dsmc.evolve(); - // Retrieve flow field data from dsmcCloud - - rhoN = dsmc.rhoN(); - rhoN.correctBoundaryConditions(); - - rhoM = dsmc.rhoM(); - rhoM.correctBoundaryConditions(); - - dsmcRhoN = dsmc.dsmcRhoN(); - dsmcRhoN.correctBoundaryConditions(); - - momentum = dsmc.momentum(); - momentum.correctBoundaryConditions(); - - linearKE = dsmc.linearKE(); - linearKE.correctBoundaryConditions(); - - internalE = dsmc.internalE(); - internalE.correctBoundaryConditions(); - - iDof = dsmc.iDof(); - iDof.correctBoundaryConditions(); - - // Retrieve surface field data from dsmcCloud - - q = dsmc.q(); - - fD = dsmc.fD(); - - // Print status of dsmcCloud - dsmc.info(); runTime.write(); diff --git a/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C b/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C index f8b831001d..9e7725ec23 100644 --- a/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C +++ b/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C @@ -44,9 +44,21 @@ int main(int argc, char *argv[]) // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + IOdictionary dsmcInitialiseDict + ( + IOobject + ( + "dsmcInitialiseDict", + mesh.time().system(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + Info<< "Initialising dsmc for Time = " << runTime.timeName() << nl << endl; - dsmcCloud dsmc("dsmc", mesh); + dsmcCloud dsmc("dsmc", mesh, dsmcInitialiseDict); label totalMolecules = dsmc.size(); diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C index 8615135a26..92a6d944d9 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C @@ -475,38 +475,99 @@ void Foam::DsmcCloud::collisions() template -void Foam::DsmcCloud::resetSurfaceDataFields() +void Foam::DsmcCloud::resetFields() { - volScalarField::GeometricBoundaryField& qBF(q_.boundaryField()); + q_ = dimensionedScalar("zero", dimensionSet(1, 0, -3, 0, 0), 0.0); - forAll(qBF, p) - { - qBF[p] = 0.0; - } - - volVectorField::GeometricBoundaryField& fDBF(fD_.boundaryField()); - - forAll(fDBF, p) - { - fDBF[p] = vector::zero; - } - - volVectorField::GeometricBoundaryField& momentumBF + fD_ = dimensionedVector ( - momentum_.boundaryField() + "zero", + dimensionSet(1, -1, -2, 0, 0), + vector::zero ); - forAll(momentumBF, p) + rhoN_ = dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL); + + rhoM_ = dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), VSMALL); + + dsmcRhoN_ = dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0); + + linearKE_ = dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0); + + internalE_ = dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0); + + iDof_ = dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL); + + momentum_ = dimensionedVector + ( + "zero", + dimensionSet(1, -2, -1, 0, 0), + vector::zero + ); +} + + +template +void Foam::DsmcCloud::calculateFields() +{ + scalarField& rhoN = rhoN_.internalField(); + rhoN = VSMALL; + + scalarField& rhoM = rhoM_.internalField(); + rhoM = VSMALL; + + scalarField& dsmcRhoN = dsmcRhoN_.internalField(); + dsmcRhoN = 0.0; + + scalarField& linearKE = linearKE_.internalField(); + linearKE = 0.0; + + scalarField& internalE = internalE_.internalField(); + internalE = 0.0; + + scalarField& iDof = iDof_.internalField(); + iDof = VSMALL; + + vectorField& momentum = momentum_.internalField(); + momentum = vector::zero; + + forAllConstIter(typename DsmcCloud, *this, iter) { - momentumBF[p] = vector::zero; + const ParcelType& p = iter(); + const label cellI = p.cell(); + + rhoN[cellI]++; + + rhoM[cellI] += constProps(p.typeId()).mass(); + + dsmcRhoN[cellI]++; + + linearKE[cellI] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U()); + + internalE[cellI] += p.Ei(); + + iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom(); + + momentum[cellI] += constProps(p.typeId()).mass()*p.U(); } - volScalarField::GeometricBoundaryField& rhoMBF(rhoM_.boundaryField()); + rhoN *= nParticle_/mesh().cellVolumes(); + rhoN_.correctBoundaryConditions(); - forAll(rhoMBF, p) - { - rhoMBF[p] = VSMALL; - } + rhoM *= nParticle_/mesh().cellVolumes(); + rhoM_.correctBoundaryConditions(); + + linearKE *= nParticle_/mesh().cellVolumes(); + linearKE_.correctBoundaryConditions(); + + internalE *= nParticle_/mesh().cellVolumes(); + internalE_.correctBoundaryConditions(); + + iDof *= nParticle_/mesh().cellVolumes(); + iDof_.correctBoundaryConditions(); + + momentum *= nParticle_/mesh().cellVolumes(); + momentum_.correctBoundaryConditions(); } @@ -542,14 +603,13 @@ template Foam::DsmcCloud::DsmcCloud ( const word& cloudName, - const volScalarField& T, - const volVectorField& U + const fvMesh& mesh ) : - Cloud(T.mesh(), cloudName, false), + Cloud(mesh, cloudName, false), DsmcBaseCloud(), cloudName_(cloudName), - mesh_(T.mesh()), + mesh_(mesh), particleProperties_ ( IOobject @@ -581,68 +641,142 @@ Foam::DsmcCloud::DsmcCloud ( IOobject ( - this->name() + "q_", + "q", mesh_.time().timeName(), mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE + IOobject::MUST_READ, + IOobject::AUTO_WRITE ), - mesh_, - dimensionedScalar("zero", dimensionSet(1, 0, -3, 0, 0), 0.0) + mesh_ ), fD_ ( IOobject ( - this->name() + "fD_", + "fD", mesh_.time().timeName(), mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE + IOobject::MUST_READ, + IOobject::AUTO_WRITE ), - mesh_, - dimensionedVector - ( - "zero", - dimensionSet(1, -1, -2, 0, 0), - vector::zero - ) + mesh_ ), - momentum_ + rhoN_ ( IOobject ( - this->name() + "momentum_", + "rhoN", mesh_.time().timeName(), mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE + IOobject::MUST_READ, + IOobject::AUTO_WRITE ), - mesh_, - dimensionedVector - ( - "zero", - dimensionSet(1, -2, -1, 0, 0), - vector::zero - ) + mesh_ ), rhoM_ ( IOobject ( - this->name() + "rhoM_", + "rhoM", mesh_.time().timeName(), mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE + IOobject::MUST_READ, + IOobject::AUTO_WRITE ), - mesh_, - dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), 0.0) + mesh_ + ), + dsmcRhoN_ + ( + IOobject + ( + "dsmcRhoN", + mesh_.time().timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh_ + ), + linearKE_ + ( + IOobject + ( + "linearKE", + mesh_.time().timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh_ + ), + internalE_ + ( + IOobject + ( + "internalE", + mesh_.time().timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh_ + ), + iDof_ + ( + IOobject + ( + "iDof", + mesh_.time().timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh_ + ), + momentum_ + ( + IOobject + ( + "momentum", + mesh_.time().timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh_ ), constProps_(), rndGen_(label(149382906) + 7183*Pstream::myProcNo()), - T_(T), - U_(U), + boundaryT_ + ( + volScalarField + ( + IOobject + ( + "boundaryT", + mesh_.time().timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh_ + ) + ), + boundaryU_ + ( + volVectorField + ( + IOobject + ( + "boundaryU", + mesh_.time().timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh_ + ) + ), binaryCollisionModel_ ( BinaryCollisionModel >::New @@ -685,7 +819,8 @@ template Foam::DsmcCloud::DsmcCloud ( const word& cloudName, - const fvMesh& mesh + const fvMesh& mesh, + const IOdictionary& dsmcInitialiseDict ) : Cloud(mesh, cloudName, false), @@ -751,6 +886,84 @@ Foam::DsmcCloud::DsmcCloud vector::zero ) ), + rhoN_ + ( + IOobject + ( + this->name() + "rhoN_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL) + ), + rhoM_ + ( + IOobject + ( + this->name() + "rhoM_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), VSMALL) + ), + dsmcRhoN_ + ( + IOobject + ( + this->name() + "dsmcRhoN_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0) + ), + linearKE_ + ( + IOobject + ( + this->name() + "linearKE_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0) + ), + internalE_ + ( + IOobject + ( + this->name() + "internalE_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0) + ), + iDof_ + ( + IOobject + ( + this->name() + "iDof_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL) + ), momentum_ ( IOobject @@ -769,28 +982,15 @@ Foam::DsmcCloud::DsmcCloud vector::zero ) ), - rhoM_ - ( - IOobject - ( - this->name() + "rhoM_", - mesh_.time().timeName(), - mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh_, - dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), 0.0) - ), constProps_(), rndGen_(label(971501) + 1526*Pstream::myProcNo()), - T_ + boundaryT_ ( volScalarField ( IOobject ( - "T", + "boundaryT", mesh_.time().timeName(), mesh_, IOobject::NO_READ, @@ -800,13 +1000,13 @@ Foam::DsmcCloud::DsmcCloud dimensionedScalar("zero", dimensionSet(0, 0, 0, 1, 0), 0.0) ) ), - U_ + boundaryU_ ( volVectorField ( IOobject ( - "U", + "boundaryU", mesh_.time().timeName(), mesh_, IOobject::NO_READ, @@ -829,18 +1029,6 @@ Foam::DsmcCloud::DsmcCloud buildConstProps(); - IOdictionary dsmcInitialiseDict - ( - IOobject - ( - "dsmcInitialiseDict", - mesh_.time().system(), - mesh_, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ); - initialise(dsmcInitialiseDict); } @@ -862,8 +1050,8 @@ void Foam::DsmcCloud::evolve() typename ParcelType::trackData td(*this); - // Reset the surface data collection fields - resetSurfaceDataFields(); + // Reset the data collection fields + resetFields(); if (debug) { @@ -878,6 +1066,9 @@ void Foam::DsmcCloud::evolve() // Calculate new velocities via stochastic collisions collisions(); + + // Calculate the volume field data + calculateFields(); } diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H index 40f3935eee..ac90e1b854 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H @@ -110,12 +110,27 @@ class DsmcCloud //- Force density at surface field volVectorField fD_; - //- Momentum density field - volVectorField momentum_; + //- number density field + volScalarField rhoN_; //- Mass density field volScalarField rhoM_; + //- dsmc particle density field + volScalarField dsmcRhoN_; + + //- linear kinetic energy density field + volScalarField linearKE_; + + //- Internal energy density field + volScalarField internalE_; + + // Internal degree of freedom density field + volScalarField iDof_; + + //- Momentum density field + volVectorField momentum_; + //- Parcel constant properties - one for each type List constProps_; @@ -127,13 +142,13 @@ class DsmcCloud scalar cachedDeltaT_; - // References to the macroscopic fields + // boundary value fields - //- Temperature - const volScalarField& T_; + //- boundary temperature + volScalarField boundaryT_; - //- Velocity - const volVectorField& U_; + //- boundary velocity + volVectorField boundaryU_; // References to the cloud sub-models @@ -165,8 +180,11 @@ class DsmcCloud //- Calculate collisions between molecules void collisions(); - //- Reset the surface data accumulation field values - void resetSurfaceDataFields(); + //- Reset the data accumulation field values to zero + void resetFields(); + + //- Calculate the volume field data + void calculateFields(); //- Disallow default bitwise copy construct DsmcCloud(const DsmcCloud&); @@ -185,21 +203,22 @@ public: // Constructors - //- Construct given name and mesh, will read Parcels from file - DsmcCloud - ( - const word& cloudName, - const volScalarField& T, - const volVectorField& U - ); - - //- Construct given name and mesh. Used to initialise. + //- Construct given name and mesh, will read Parcels and fields from + // file DsmcCloud ( const word& cloudName, const fvMesh& mesh ); + //- Construct given name, mesh and initialisation dictionary. + DsmcCloud + ( + const word& cloudName, + const fvMesh& mesh, + const IOdictionary& dsmcInitialiseDict + ); + //- Destructor virtual ~DsmcCloud(); @@ -268,20 +287,48 @@ public: //- Return non-const force density at boundary field reference inline volVectorField::GeometricBoundaryField& fDBF(); - //- Return non-const momentum density boundary field reference - inline volVectorField::GeometricBoundaryField& momentumBF(); - //- Return non-const mass density boundary field reference inline volScalarField::GeometricBoundaryField& rhoMBF(); + //- Return non-const momentum density boundary field reference + inline volVectorField::GeometricBoundaryField& momentumBF(); + // References to the macroscopic fields //- Return macroscopic temperature - inline const volScalarField& T() const; + inline const volScalarField& boundaryT() const; //- Return macroscopic velocity - inline const volVectorField& U() const; + inline const volVectorField& boundaryU() const; + + //- Return heat flux at surface field + inline const volScalarField& q() const; + + //- Return force density at surface field + inline const volVectorField& fD() const; + + //- Return the real particle number density field + inline const volScalarField& rhoN() const; + + //- Return the particle mass density field + inline const volScalarField& rhoM() const; + + //- Return the field of number of DSMC particles + inline const volScalarField& dsmcRhoN() const; + + //- Return the total linear kinetic energy (translational and + // thermal density field + inline const volScalarField& linearKE() const; + + //- Return the internal energy density field + inline const volScalarField& internalE() const; + + //- Return the average internal degrees of freedom field + inline const volScalarField& iDof() const; + + //- Return the momentum density field + inline const volVectorField& momentum() const; // Kinetic theory helper functions @@ -395,35 +442,6 @@ public: void dumpParticlePositions() const; - // Fields - - //- Return heat flux at surface field - inline const volScalarField& q() const; - - //- Return force density at surface field - inline const volVectorField& fD() const; - - //- Return the real particle number density field - inline const tmp rhoN() const; - - //- Return the particle mass density field - inline const volScalarField& rhoM(); - - //- Return the momentum density field - inline const volVectorField& momentum(); - - //- 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; - - //- Return the field of number of DSMC particles - inline const tmp dsmcRhoN() const; // Cloud evolution functions @@ -441,7 +459,6 @@ public: //- Evolve the cloud (move, collide) void evolve(); - //- Clear the Cloud inline void clear(); diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H index 07ac87bd2c..5b780ee171 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H @@ -167,16 +167,18 @@ Foam::DsmcCloud::rhoMBF() template -inline const Foam::volScalarField& Foam::DsmcCloud::T() const +inline const Foam::volScalarField& +Foam::DsmcCloud::boundaryT() const { - return T_; + return boundaryT_; } template -inline const Foam::volVectorField& Foam::DsmcCloud::U() const +inline const Foam::volVectorField& +Foam::DsmcCloud::boundaryU() const { - return U_; + return boundaryU_; } @@ -392,228 +394,56 @@ inline const Foam::volVectorField& Foam::DsmcCloud::fD() const template -inline const Foam::tmp +inline const Foam::volScalarField& Foam::DsmcCloud::rhoN() const { - tmp trhoN - ( - new volScalarField - ( - IOobject - ( - this->name() + "rhoN", - this->db().time().timeName(), - this->db(), - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh_, - dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL) - ) - ); - - 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; + return rhoN_; } template -inline const Foam::volScalarField& Foam::DsmcCloud::rhoM() +inline const Foam::volScalarField& Foam::DsmcCloud::rhoM() const { - scalarField& rM = rhoM_.internalField(); - - rM = VSMALL; - - forAllConstIter(typename DsmcCloud, *this, iter) - { - const ParcelType& p = iter(); - const label cellI = p.cell(); - - rM[cellI] += constProps(p.typeId()).mass(); - } - - rM *= nParticle_/mesh().cellVolumes(); - return rhoM_; } template -inline const Foam::volVectorField& Foam::DsmcCloud::momentum() -{ - vectorField& mom = momentum_.internalField(); - - mom = vector::zero; - - forAllConstIter(typename DsmcCloud, *this, iter) - { - const ParcelType& p = iter(); - const label cellI = p.cell(); - - mom[cellI] += constProps(p.typeId()).mass()*p.U(); - } - - mom *= nParticle_/mesh().cellVolumes(); - - return momentum_; -} - - -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", dimensionSet(0, -3, 0, 0, 0), VSMALL) - ) - ); - - scalarField& iDof = tiDof().internalField(); - - forAllConstIter(typename DsmcCloud, *this, iter) - { - const ParcelType& p = iter(); - const label cellI = p.cell(); - - iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom(); - } - - iDof *= nParticle_/mesh().cellVolumes(); - - return tiDof; -} - - -template -inline const Foam::tmp +inline const Foam::volScalarField& Foam::DsmcCloud::dsmcRhoN() const { - tmp tdsmcRhoN - ( - new volScalarField - ( - IOobject - ( - this->name() + "dsmcRhoN", - this->db().time().timeName(), - this->db(), - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh_, - dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0) - ) - ); + return dsmcRhoN_; +} - scalarField& dsmcRhoN = tdsmcRhoN().internalField(); - forAllConstIter(typename DsmcCloud, *this, iter) - { - const ParcelType& p = iter(); - const label cellI = p.cell(); - dsmcRhoN[cellI]++; - } +template +inline const Foam::volScalarField& +Foam::DsmcCloud::linearKE() const +{ + return linearKE_; +} - return tdsmcRhoN; + +template +inline const Foam::volScalarField& +Foam::DsmcCloud::internalE() const +{ + return internalE_; +} + + +template +inline const Foam::volScalarField& +Foam::DsmcCloud::iDof() const +{ + return iDof_; +} + + +template +inline const Foam::volVectorField& Foam::DsmcCloud::momentum() const +{ + return momentum_; } diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C index 425e47a705..1f45f25e49 100644 --- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C +++ b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C @@ -39,11 +39,10 @@ namespace Foam Foam::dsmcCloud::dsmcCloud ( const word& cloudName, - const volScalarField& T, - const volVectorField& U + const fvMesh& mesh ) : - DsmcCloud(cloudName, T, U) + DsmcCloud(cloudName, mesh) { dsmcParcel::readFields(*this); } @@ -52,10 +51,11 @@ Foam::dsmcCloud::dsmcCloud Foam::dsmcCloud::dsmcCloud ( const word& cloudName, - const fvMesh& mesh + const fvMesh& mesh, + const IOdictionary& dsmcInitialiseDict ) : - DsmcCloud(cloudName, mesh) + DsmcCloud(cloudName, mesh, dsmcInitialiseDict) {} diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H index ad63aca9fc..953e487e8d 100644 --- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H +++ b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H @@ -69,21 +69,22 @@ public: // Constructors - //- Construct from components - dsmcCloud - ( - const word& cloudName, - const volScalarField& T, - const volVectorField& U - ); - - //- Construct from name and mesh, used to initialise. + //- Construct given name and mesh, will read Parcels and fields from + // file dsmcCloud ( const word& cloudName, const fvMesh& mesh ); + //- Construct given name, mesh and initialisation dictionary. + dsmcCloud + ( + const word& cloudName, + const fvMesh& mesh, + const IOdictionary& dsmcInitialiseDict + ); + //- Destructor ~dsmcCloud(); diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C index 0f0e135e4c..f02fcad0dc 100644 --- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C +++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C @@ -136,12 +136,12 @@ void Foam::FreeStream::inflow() const volScalarField::GeometricBoundaryField& boundaryT ( - cloud.T().boundaryField() + cloud.boundaryT().boundaryField() ); const volVectorField::GeometricBoundaryField& boundaryU ( - cloud.U().boundaryField() + cloud.boundaryU().boundaryField() ); forAll(patches_, p) @@ -165,7 +165,8 @@ void Foam::FreeStream::inflow() if (min(boundaryT[patchI]) < SMALL) { FatalErrorIn ("Foam::FreeStream::inflow()") - << "Zero boundary temperature detected, check boundaryT condition." << nl + << "Zero boundary temperature detected, check boundaryT " + << "condition." << nl << nl << abort(FatalError); } diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C index 6ad452793b..51452e7c1f 100644 --- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C +++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C @@ -101,7 +101,7 @@ void Foam::MaxwellianThermal::correct // Other tangential unit vector vector tw2 = nw^tw1; - scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace]; + scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace]; scalar mass = cloud.constProps(typeId).mass(); @@ -115,7 +115,7 @@ void Foam::MaxwellianThermal::correct - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw ); - U += cloud.U().boundaryField()[wppIndex][wppLocalFace]; + U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace]; Ei = cloud.equipartitionInternalEnergy(T, iDof); } diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C index d491dfa406..9d22415f3c 100644 --- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C +++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C @@ -106,7 +106,7 @@ void Foam::MixedDiffuseSpecular::correct // Other tangential unit vector vector tw2 = nw^tw1; - scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace]; + scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace]; scalar mass = cloud.constProps(typeId).mass(); @@ -120,7 +120,7 @@ void Foam::MixedDiffuseSpecular::correct - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw ); - U += cloud.U().boundaryField()[wppIndex][wppLocalFace]; + U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace]; Ei = cloud.equipartitionInternalEnergy(T, iDof); } From c0972656454cd443c5c26a689204ed049f1ecdad Mon Sep 17 00:00:00 2001 From: graham Date: Tue, 21 Jul 2009 15:10:23 +0100 Subject: [PATCH 05/16] Adding normal velocity weighted measurements to face values of linearKE, internalE, rhoN and iDof fields to measure surface temperatures. --- .../clouds/Templates/DsmcCloud/DsmcCloud.C | 7 ---- .../clouds/Templates/DsmcCloud/DsmcCloud.H | 15 ++++++++ .../clouds/Templates/DsmcCloud/DsmcCloudI.H | 38 +++++++++++++++++-- .../parcels/Templates/DsmcParcel/DsmcParcel.C | 38 +++++++++++++++---- 4 files changed, 80 insertions(+), 18 deletions(-) diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C index 92a6d944d9..bb51dbd24b 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C @@ -511,25 +511,18 @@ template void Foam::DsmcCloud::calculateFields() { scalarField& rhoN = rhoN_.internalField(); - rhoN = VSMALL; scalarField& rhoM = rhoM_.internalField(); - rhoM = VSMALL; scalarField& dsmcRhoN = dsmcRhoN_.internalField(); - dsmcRhoN = 0.0; scalarField& linearKE = linearKE_.internalField(); - linearKE = 0.0; scalarField& internalE = internalE_.internalField(); - internalE = 0.0; scalarField& iDof = iDof_.internalField(); - iDof = VSMALL; vectorField& momentum = momentum_.internalField(); - momentum = vector::zero; forAllConstIter(typename DsmcCloud, *this, iter) { diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H index ac90e1b854..67deff457b 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H @@ -287,9 +287,24 @@ public: //- Return non-const force density at boundary field reference inline volVectorField::GeometricBoundaryField& fDBF(); + //- Return non-const number density boundary field reference + inline volScalarField::GeometricBoundaryField& rhoNBF(); + //- Return non-const mass density boundary field reference inline volScalarField::GeometricBoundaryField& rhoMBF(); + //- Return non-const linear kinetic energy density boundary + // field reference + inline volScalarField::GeometricBoundaryField& linearKEBF(); + + //- Return non-const internal energy density boundary field + // reference + inline volScalarField::GeometricBoundaryField& internalEBF(); + + //- Return non-const internal degree of freedom density boundary + // field reference + inline volScalarField::GeometricBoundaryField& iDofBF(); + //- Return non-const momentum density boundary field reference inline volVectorField::GeometricBoundaryField& momentumBF(); diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H index 5b780ee171..f2dcb16bce 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H @@ -151,10 +151,10 @@ Foam::DsmcCloud::fDBF() template -inline Foam::volVectorField::GeometricBoundaryField& -Foam::DsmcCloud::momentumBF() +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::rhoNBF() { - return momentum_.boundaryField(); + return rhoN_.boundaryField(); } @@ -166,6 +166,38 @@ Foam::DsmcCloud::rhoMBF() } +template +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::linearKEBF() +{ + return linearKE_.boundaryField(); +} + + +template +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::internalEBF() +{ + return internalE_.boundaryField(); +} + + +template +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::iDofBF() +{ + return iDof_.boundaryField(); +} + + +template +inline Foam::volVectorField::GeometricBoundaryField& +Foam::DsmcCloud::momentumBF() +{ + return momentum_.boundaryField(); +} + + template inline const Foam::volScalarField& Foam::DsmcCloud::boundaryT() const diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C index e11a5e73ea..ffd3485fdb 100644 --- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C +++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C @@ -130,13 +130,24 @@ void Foam::DsmcParcel::hitWallPatch nw /= mag(nw); scalar U_dot_nw = U_ & nw; + vector Ut = U_ - U_dot_nw*nw; - td.cloud().momentumBF()[wppIndex][wppLocalFace] += - m*Ut/max(mag(U_dot_nw)*fA, VSMALL); + scalar invMagUnfA = 1/max(mag(U_dot_nw)*fA, VSMALL); - td.cloud().rhoMBF()[wppIndex][wppLocalFace] += - m/max(mag(U_dot_nw)*fA, VSMALL); + td.cloud().rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA; + + td.cloud().rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA; + + td.cloud().linearKEBF()[wppIndex][wppLocalFace] += + 0.5*m*(U_ & U_)*invMagUnfA; + + td.cloud().internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA; + + td.cloud().iDofBF()[wppIndex][wppLocalFace] += + constProps.internalDegreesOfFreedom()*invMagUnfA; + + td.cloud().momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA; // pre-interaction energy scalar preIE = 0.5*m*(U_ & U_) + Ei_; @@ -154,13 +165,24 @@ void Foam::DsmcParcel::hitWallPatch ); U_dot_nw = U_ & nw; + Ut = U_ - U_dot_nw*nw; - td.cloud().momentumBF()[wppIndex][wppLocalFace] += - m*Ut/max(mag(U_dot_nw)*fA, VSMALL); + invMagUnfA = 1/max(mag(U_dot_nw)*fA, VSMALL); - td.cloud().rhoMBF()[wppIndex][wppLocalFace] += - m/max(mag(U_dot_nw)*fA, VSMALL); + td.cloud().rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA; + + td.cloud().rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA; + + td.cloud().linearKEBF()[wppIndex][wppLocalFace] += + 0.5*m*(U_ & U_)*invMagUnfA; + + td.cloud().internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA; + + td.cloud().iDofBF()[wppIndex][wppLocalFace] += + constProps.internalDegreesOfFreedom()*invMagUnfA; + + td.cloud().momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA; // post-interaction energy scalar postIE = 0.5*m*(U_ & U_) + Ei_; From 85aa1c889047e5443da1d5b3cdadf8ad81036525 Mon Sep 17 00:00:00 2001 From: graham Date: Tue, 21 Jul 2009 17:17:01 +0100 Subject: [PATCH 06/16] Reordering of internal temperature calculation to avoid FPEs. iDof was defaulting to VSMALL (1e-300) for monoatomics but being multiplied by kb: O(1e-23) for field calculation. --- .../functionObjects/utilities/dsmcFields/dsmcFields.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C index 7056fba521..ab6190af8b 100644 --- a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C +++ b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C @@ -184,7 +184,7 @@ void Foam::dsmcFields::write() obr_, IOobject::NO_READ ), - 2.0/(dsmcCloud::kb*iDofMean)*internalEMean + (2.0/dsmcCloud::kb)*(internalEMean/iDofMean) ); Info<< " Calculating overallT field." << endl; From 8e885770126d4eef80a8ecf37e75bb8af86cedc2 Mon Sep 17 00:00:00 2001 From: graham Date: Wed, 19 Aug 2009 15:20:00 +0100 Subject: [PATCH 07/16] Revamping distribution class to use for dsmc tests. --- etc/controlDict | 1 + .../distribution/distribution.C | 108 ++++++++++++------ .../distribution/distribution.H | 21 +++- 3 files changed, 90 insertions(+), 40 deletions(-) diff --git a/etc/controlDict b/etc/controlDict index c1d7416515..6f47c9c273 100644 --- a/etc/controlDict +++ b/etc/controlDict @@ -377,6 +377,7 @@ DebugSwitches displacementLaplacian 0; displacementSBRStress 0; distanceSurface 0; + distribution 0; downwind 0; dragModel 0; duplicatePoints 0; diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C index 439ed64bdf..7cc9e70c0e 100644 --- a/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C +++ b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C @@ -25,27 +25,49 @@ License \*----------------------------------------------------------------------------*/ #include "distribution.H" +#include "OFstream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { + defineTypeNameAndDebug(distribution, 0); +} + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +void Foam::distribution::write +( + const fileName& file, + const List >& pairs +) +{ + OFstream os(file); + + forAll(pairs, i) + { + os << pairs[i].first() << ' ' << pairs[i].second() << nl; + } +} + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -distribution::distribution() +Foam::distribution::distribution() : Map