From f6e6c62babff362b0b2f0936f5e7884ed85e4f07 Mon Sep 17 00:00:00 2001 From: henry Date: Wed, 11 Nov 2009 15:22:11 +0000 Subject: [PATCH] Added measurement of surface values of velocity and temperature to allow slip and jump measurement. Added measurement of volume and surface pressure. Cases may now also be run as 1D or 2D by using empty patches. --- .../discreteMethods/dsmc/dsmcFoam/createFields.H | 162 --------- .../discreteMethods/dsmc/dsmcFoam/dsmcFoam.C | 42 +-- .../miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C | 11 +- .../preProcessing/dsmcInitialise/dsmcInitialise.C | 14 +- src/lagrangian/dsmc/Make/options | 6 +- .../dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C | 365 ++++++++++++++++++--- .../dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H | 148 +++++---- .../dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H | 309 ++++------------- .../dsmc/clouds/derived/dsmcCloud/dsmcCloud.C | 16 +- .../dsmc/clouds/derived/dsmcCloud/dsmcCloud.H | 12 +- .../dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C | 86 ++++- .../dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H | 1 - .../parcels/Templates/DsmcParcel/DsmcParcelI.H | 2 + .../parcels/Templates/DsmcParcel/DsmcParcelIO.C | 10 +- .../makeDsmcParcelWallInteractionModels.C | 7 + .../InflowBoundaryModel/FreeStream/FreeStream.C | 9 +- .../MaxwellianThermal/MaxwellianThermal.C | 4 +- .../MixedDiffuseSpecular.C} | 92 +++--- .../MixedDiffuseSpecular/MixedDiffuseSpecular.H} | 73 +++-- .../utilities/dsmcFields/dsmcFields.C | 43 ++- 20 files changed, 745 insertions(+), 667 deletions(-) delete mode 100644 applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H copy src/lagrangian/dsmc/submodels/WallInteractionModel/{MaxwellianThermal/MaxwellianThermal.C => MixedDiffuseSpecular/MixedDiffuseSpecular.C} (52%) copy src/lagrangian/dsmc/{clouds/derived/dsmcCloud/dsmcCloud.H => submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H} (62%) diff --git a/applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H b/applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H deleted file mode 100644 index d024bd20..00000000 --- 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 84b49945..f5a57f15 100644 --- a/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C +++ b/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C @@ -41,53 +41,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/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C b/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C index 95077df2..34f23185 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/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C b/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C index f8b83100..9e7725ec 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/Make/options b/src/lagrangian/dsmc/Make/options index 15874b7b..55865dfa 100644 --- a/src/lagrangian/dsmc/Make/options +++ b/src/lagrangian/dsmc/Make/options @@ -1,7 +1,9 @@ EXE_INC = \ -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/lagrangian/basic/lnInclude + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude LIB_LIBS = \ -llagrangian \ - -lfiniteVolume + -lfiniteVolume \ + -lmeshTools diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C index 23aa3dc4..a3f7093c 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C @@ -299,7 +299,7 @@ void Foam::DsmcCloud::collisions() // Temporary storage for subCells List > subCells(8); - scalar deltaT = cachedDeltaT(); + scalar deltaT = mesh().time().deltaTValue(); label collisionCandidates = 0; @@ -475,21 +475,92 @@ 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); + + fD_ = dimensionedVector + ( + "zero", + dimensionSet(1, -1, -2, 0, 0), + vector::zero + ); + + 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 + ); +} - forAll(qBF, p) - { - qBF[p] = 0.0; - } - volVectorField::GeometricBoundaryField& fDBF(fD_.boundaryField()); +template +void Foam::DsmcCloud::calculateFields() +{ + scalarField& rhoN = rhoN_.internalField(); + + scalarField& rhoM = rhoM_.internalField(); + + scalarField& dsmcRhoN = dsmcRhoN_.internalField(); + + scalarField& linearKE = linearKE_.internalField(); + + scalarField& internalE = internalE_.internalField(); + + scalarField& iDof = iDof_.internalField(); + + vectorField& momentum = momentum_.internalField(); - forAll(fDBF, p) + forAllConstIter(typename DsmcCloud, *this, iter) { - fDBF[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(); } + + rhoN *= nParticle_/mesh().cellVolumes(); + rhoN_.correctBoundaryConditions(); + + 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(); } @@ -525,14 +596,14 @@ template Foam::DsmcCloud::DsmcCloud ( const word& cloudName, - const volScalarField& T, - const volVectorField& U + const fvMesh& mesh, + bool readFields ) : - Cloud(T.mesh(), cloudName, false), + Cloud(mesh, cloudName, false), DsmcBaseCloud(), cloudName_(cloudName), - mesh_(T.mesh()), + mesh_(mesh), particleProperties_ ( IOobject @@ -564,37 +635,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 + mesh_ + ), + rhoN_ + ( + IOobject ( - "zero", - dimensionSet(1, -1, -2, 0, 0), - vector::zero - ) + "rhoN", + mesh_.time().timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh_ + ), + rhoM_ + ( + IOobject + ( + "rhoM", + mesh_.time().timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + 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 @@ -637,7 +813,8 @@ template Foam::DsmcCloud::DsmcCloud ( const word& cloudName, - const fvMesh& mesh + const fvMesh& mesh, + const IOdictionary& dsmcInitialiseDict ) : Cloud(mesh, cloudName, false), @@ -703,15 +880,111 @@ 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 + ( + this->name() + "momentum_", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedVector + ( + "zero", + dimensionSet(1, -2, -1, 0, 0), + vector::zero + ) + ), constProps_(), rndGen_(label(971501) + 1526*Pstream::myProcNo()), - T_ + boundaryT_ ( volScalarField ( IOobject ( - "T", + "boundaryT", mesh_.time().timeName(), mesh_, IOobject::NO_READ, @@ -721,13 +994,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, @@ -750,18 +1023,6 @@ Foam::DsmcCloud::DsmcCloud buildConstProps(); - IOdictionary dsmcInitialiseDict - ( - IOobject - ( - "dsmcInitialiseDict", - mesh_.time().system(), - mesh_, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ); - initialise(dsmcInitialiseDict); } @@ -778,13 +1039,10 @@ Foam::DsmcCloud::~DsmcCloud() template void Foam::DsmcCloud::evolve() { - // cache the value of deltaT for this timestep - storeDeltaT(); - typename ParcelType::trackData td(*this); - // Reset the surface data collection fields - resetSurfaceDataFields(); + // Reset the data collection fields + resetFields(); if (debug) { @@ -799,6 +1057,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 76b1dd7e..2a9088d5 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H @@ -110,24 +110,41 @@ class DsmcCloud //- Force density at surface field volVectorField fD_; + //- 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_; //- Random number generator Random rndGen_; - //- In-cloud cache of deltaT, lookup in submodels and parcel is - // expensive - scalar cachedDeltaT_; + // boundary value fields - // References to the macroscopic fields + //- boundary temperature + volScalarField boundaryT_; - //- Temperature - const volScalarField& T_; - - //- Velocity - const volVectorField& U_; + //- boundary velocity + volVectorField boundaryU_; // References to the cloud sub-models @@ -159,8 +176,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&); @@ -179,19 +199,21 @@ public: // Constructors - //- Construct given name and mesh, will read Parcels from file + //- Construct given name and mesh, will read Parcels and fields from + // file DsmcCloud ( const word& cloudName, - const volScalarField& T, - const volVectorField& U + const fvMesh& mesh, + bool readFields = true ); - //- Construct given name and mesh. Used to initialise. + //- Construct given name, mesh and initialisation dictionary. DsmcCloud ( const word& cloudName, - const fvMesh& mesh + const fvMesh& mesh, + const IOdictionary& dsmcInitialiseDict ); @@ -247,35 +269,72 @@ public: //- Return refernce to the random object inline Random& rndGen(); - //- Store (cache) the current value of deltaT - inline void storeDeltaT(); - //- Return the cached value of deltaT - inline scalar cachedDeltaT() const; + // References to the boundary fields for surface data collection + //- Return non-const heat flux boundary field reference + inline volScalarField::GeometricBoundaryField& qBF(); - // References to the surface data collection fields + //- Return non-const force density at boundary field reference + inline volVectorField::GeometricBoundaryField& fDBF(); - //- Return heat flux at surface field - inline const volScalarField& q() const; + //- Return non-const number density boundary field reference + inline volScalarField::GeometricBoundaryField& rhoNBF(); - //- Return non-const heat flux at surface field - inline volScalarField& q(); + //- Return non-const mass density boundary field reference + inline volScalarField::GeometricBoundaryField& rhoMBF(); - //- Return force density at surface field - inline const volVectorField& fD() const; + //- 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 force density at surface field - inline volVectorField& fD(); + //- 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(); // 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 @@ -338,6 +397,7 @@ public: scalar mass ) const; + // Sub-models //- Return reference to binary elastic collision model @@ -389,29 +449,6 @@ public: void dumpParticlePositions() const; - // Fields - - //- 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; - - //- 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 @@ -429,11 +466,8 @@ 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 3529702f..afb7f642 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H @@ -121,58 +121,82 @@ inline Foam::Random& Foam::DsmcCloud::rndGen() template -inline void Foam::DsmcCloud::storeDeltaT() +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::qBF() { - cachedDeltaT_ = mesh().time().deltaT().value(); + return q_.boundaryField(); } template -inline Foam::scalar Foam::DsmcCloud::cachedDeltaT() const +inline Foam::volVectorField::GeometricBoundaryField& +Foam::DsmcCloud::fDBF() { - return cachedDeltaT_; + return fD_.boundaryField(); } template -inline const Foam::volScalarField& Foam::DsmcCloud::q() const +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::rhoNBF() { - return q_; + return rhoN_.boundaryField(); } template -inline Foam::volScalarField& Foam::DsmcCloud::q() +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::rhoMBF() { - return q_; + return rhoM_.boundaryField(); } template -inline const Foam::volVectorField& Foam::DsmcCloud::fD() const +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::linearKEBF() { - return fD_; + return linearKE_.boundaryField(); } template -inline Foam::volVectorField& Foam::DsmcCloud::fD() +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::internalEBF() { - return fD_; + return internalE_.boundaryField(); +} + + +template +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::iDofBF() +{ + return iDof_.boundaryField(); } template -inline const Foam::volScalarField& Foam::DsmcCloud::T() const +inline Foam::volVectorField::GeometricBoundaryField& +Foam::DsmcCloud::momentumBF() { - return T_; + return momentum_.boundaryField(); } template -inline const Foam::volVectorField& Foam::DsmcCloud::U() const +inline const Foam::volScalarField& +Foam::DsmcCloud::boundaryT() const { - return U_; + return boundaryT_; +} + + +template +inline const Foam::volVectorField& +Foam::DsmcCloud::boundaryU() const +{ + return boundaryU_; } @@ -312,7 +336,8 @@ inline Foam::scalar Foam::DsmcCloud::maxwellianAverageSpeed scalar mass ) const { - return 2.0*sqrt(2.0*kb*temperature/(mathematicalConstant::pi*mass)); + return + 2.0*sqrt(2.0*kb*temperature/(mathematicalConstant::pi*mass)); } @@ -323,7 +348,8 @@ inline Foam::scalarField Foam::DsmcCloud::maxwellianAverageSpeed scalar mass ) const { - return 2.0*sqrt(2.0*kb*temperature/(mathematicalConstant::pi*mass)); + return + 2.0*sqrt(2.0*kb*temperature/(mathematicalConstant::pi*mass)); } @@ -374,265 +400,70 @@ Foam::DsmcCloud::maxwellianMostProbableSpeed template -inline const Foam::tmp -Foam::DsmcCloud::rhoN() const +inline const Foam::volScalarField& Foam::DsmcCloud::q() 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 q_; } template -inline const Foam::tmp -Foam::DsmcCloud::rhoM() const +inline const Foam::volVectorField& Foam::DsmcCloud::fD() const { - 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& 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; + return fD_; } template -inline const Foam::tmp -Foam::DsmcCloud::dsmcRhoN() const +inline const Foam::volScalarField& +Foam::DsmcCloud::rhoN() 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; + return rhoN_; } template -inline const Foam::tmp -Foam::DsmcCloud::momentum() const +inline const Foam::volScalarField& Foam::DsmcCloud::rhoM() 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(); - } + return rhoM_; +} - momentum *= nParticle_/mesh().cellVolumes(); - return tmomentum; +template +inline const Foam::volScalarField& +Foam::DsmcCloud::dsmcRhoN() const +{ + return dsmcRhoN_; } template -inline const Foam::tmp +inline const Foam::volScalarField& 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; + return linearKE_; } template -inline const Foam::tmp +inline const Foam::volScalarField& 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; + return internalE_; } template -inline const Foam::tmp +inline const Foam::volScalarField& 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(); - } + return iDof_; +} - iDof *= nParticle_/mesh().cellVolumes(); - return tiDof; +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 425e47a7..9d1510c5 100644 --- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C +++ b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C @@ -39,23 +39,27 @@ namespace Foam Foam::dsmcCloud::dsmcCloud ( const word& cloudName, - const volScalarField& T, - const volVectorField& U + const fvMesh& mesh, + bool readFields ) : - DsmcCloud(cloudName, T, U) + DsmcCloud(cloudName, mesh, readFields) { - dsmcParcel::readFields(*this); + if (readFields) + { + dsmcParcel::readFields(*this); + } } 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 ad63aca9..07280e16 100644 --- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H +++ b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H @@ -69,19 +69,21 @@ public: // Constructors - //- Construct from components + //- Construct given name and mesh, will read Parcels and fields from + // file dsmcCloud ( const word& cloudName, - const volScalarField& T, - const volVectorField& U + const fvMesh& mesh, + bool readFields = true ); - //- Construct from name and mesh, used to initialise. + //- Construct given name, mesh and initialisation dictionary. dsmcCloud ( const word& cloudName, - const fvMesh& mesh + const fvMesh& mesh, + const IOdictionary& dsmcInitialiseDict ); //- Destructor diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C index 3d409b5c..3d752aa9 100644 --- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C +++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C @@ -25,7 +25,7 @@ License \*---------------------------------------------------------------------------*/ #include "DsmcParcel.H" -#include "dimensionedConstants.H" +#include "meshTools.H" // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // @@ -44,16 +44,31 @@ bool Foam::DsmcParcel::move const polyMesh& mesh = td.cloud().pMesh(); const polyBoundaryMesh& pbMesh = mesh.boundaryMesh(); - const scalar deltaT = td.cloud().cachedDeltaT(); + const scalar deltaT = mesh.time().deltaTValue(); scalar tEnd = (1.0 - p.stepFraction())*deltaT; const scalar dtMax = tEnd; + // For reduced-D cases, the velocity used to track needs to be + // constrained, but the actual U_ of the parcel must not be + // altered or used, as it is altered by patch interactions an + // needs to retain its 3D value for collision purposes. + vector Utracking = U_; + while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL) { + // Apply correction to position for reduced-D cases + meshTools::constrainToMeshCentre(mesh, p.position()); + + Utracking = U_; + + // Apply correction to velocity to constrain tracking for + // reduced-D cases + meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking); + // Set the Lagrangian time-step scalar dt = min(dtMax, tEnd); - dt *= p.trackToFace(p.position() + dt*U_, td); + dt *= p.trackToFace(p.position() + dt*Utracking, td); tEnd -= dt; @@ -114,10 +129,41 @@ 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().pMesh().time().deltaTValue(); + 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; + + scalar invMagUnfA = 1/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_; @@ -133,27 +179,40 @@ void Foam::DsmcParcel::hitWallPatch typeId_ ); - // post-interaction energy - scalar postIE = 0.5*m*(U_ & U_) + Ei_; + U_dot_nw = U_ & nw; - // post-interaction momentum - vector postIMom = m*U_; + Ut = U_ - U_dot_nw*nw; - label wppIndex = wpp.index(); + invMagUnfA = 1/max(mag(U_dot_nw)*fA, VSMALL); - label wppLocalFace = wpp.whichFace(this->face()); + td.cloud().rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA; - const scalar fA = mag(wpp.faceAreas()[wppLocalFace]); + 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; - const scalar deltaT = td.cloud().cachedDeltaT(); + 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_; + + // post-interaction momentum + vector postIMom = m*U_; 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; } @@ -213,4 +272,3 @@ void Foam::DsmcParcel::transformProperties #include "DsmcParcelIO.C" // ************************************************************************* // - diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H index af8e350e..d96d5158 100644 --- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H +++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.H @@ -42,7 +42,6 @@ SourceFiles #include "IOstream.H" #include "autoPtr.H" #include "contiguous.H" -#include "mathematicalConstants.H" #include "DsmcCloud.H" diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H index 239571ce..372a3537 100644 --- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H +++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelI.H @@ -24,6 +24,8 @@ License \*---------------------------------------------------------------------------*/ +#include "mathematicalConstants.H" + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C index 66775910..598c01c9 100644 --- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C +++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcelIO.C @@ -74,10 +74,7 @@ Foam::DsmcParcel::DsmcParcel template -void Foam::DsmcParcel::readFields -( - DsmcCloud& c -) +void Foam::DsmcParcel::readFields(DsmcCloud& c) { if (!c.size()) { @@ -107,10 +104,7 @@ void Foam::DsmcParcel::readFields template -void Foam::DsmcParcel::writeFields -( - const DsmcCloud& c -) +void Foam::DsmcParcel::writeFields(const DsmcCloud& c) { Particle::writeFields(c); diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C index f15d0bda..b67f66ec 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/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C index 17167b50..5fa77347 100644 --- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C +++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C @@ -126,7 +126,7 @@ void Foam::FreeStream::inflow() const polyMesh& mesh(cloud.mesh()); - const scalar deltaT = cloud.cachedDeltaT(); + const scalar deltaT = mesh.time().deltaTValue(); Random& rndGen(cloud.rndGen()); @@ -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 e7213561..f29ed04e 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/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C similarity index 52% copy from src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C copy to src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C index e7213561..e2484140 100644 --- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C +++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C @@ -24,32 +24,33 @@ License \*---------------------------------------------------------------------------*/ -#include "MaxwellianThermal.H" +#include "MixedDiffuseSpecular.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -Foam::MaxwellianThermal::MaxwellianThermal +Foam::MixedDiffuseSpecular::MixedDiffuseSpecular ( const dictionary& dict, CloudType& cloud ) : - WallInteractionModel(dict, cloud, typeName) + WallInteractionModel(dict, cloud, typeName), + diffuseFraction_(readScalar(this->coeffDict().lookup("diffuseFraction"))) {} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template -Foam::MaxwellianThermal::~MaxwellianThermal() +Foam::MixedDiffuseSpecular::~MixedDiffuseSpecular() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -void Foam::MaxwellianThermal::correct +void Foam::MixedDiffuseSpecular::correct ( const wallPolyPatch& wpp, const label faceId, @@ -68,56 +69,71 @@ void Foam::MaxwellianThermal::correct nw /= mag(nw); // Normal velocity magnitude - scalar magUn = U & nw; - - // Wall tangential velocity (flow direction) - vector Ut = U - magUn*nw; + scalar U_dot_nw = U & nw; CloudType& cloud(this->owner()); Random& rndGen(cloud.rndGen()); - while (mag(Ut) < SMALL) + if (diffuseFraction_ > rndGen.scalar01()) { - // 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. + // Diffuse reflection - 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()) - ); + // Wall tangential velocity (flow direction) + vector Ut = U - U_dot_nw*nw; - magUn = U & 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. - Ut = U - magUn*nw; - } + 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()) + ); + + U_dot_nw = U & nw; + + Ut = U - U_dot_nw*nw; + } + + // Wall tangential unit vector + vector tw1 = Ut/mag(Ut); - // Wall tangential unit vector - vector tw1 = Ut/mag(Ut); + // Other tangential unit vector + vector tw2 = nw^tw1; - // Other tangential unit vector - vector tw2 = nw^tw1; + scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace]; - scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace]; + scalar mass = cloud.constProps(typeId).mass(); - scalar mass = cloud.constProps(typeId).mass(); + scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom(); - 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 = - sqrt(CloudType::kb*T/mass) - *( - rndGen.GaussNormal()*tw1 - + rndGen.GaussNormal()*tw2 - - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw - ); + U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace]; - U += cloud.U().boundaryField()[wppIndex][wppLocalFace]; + Ei = cloud.equipartitionInternalEnergy(T, iDof); + } + else + { + // Specular reflection + + if (U_dot_nw > 0.0) + { + U -= 2.0*U_dot_nw*nw; + } + } - Ei = cloud.equipartitionInternalEnergy(T, iDof); } diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H similarity index 62% copy from src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H copy to src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H index ad63aca9..b84f43d8 100644 --- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H +++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H @@ -23,83 +23,84 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - Foam::dsmcCloud + Foam::MixedDiffuseSpecular Description - Cloud class to simulate dsmc parcels - -SourceFiles - dsmcCloud.C + Wall interaction setting microscopic velocity to a random one drawn from a + Maxwellian distribution corresponding to a specified temperature \*---------------------------------------------------------------------------*/ -#ifndef dsmcCloud_H -#define dsmcCloud_H +#ifndef MixedDiffuseSpecular_H +#define MixedDiffuseSpecular_H -#include "DsmcCloud.H" -#include "dsmcParcel.H" +#include "WallInteractionModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { - /*---------------------------------------------------------------------------*\ - Class dsmcCloud Declaration + Class MixedDiffuseSpecular Declaration \*---------------------------------------------------------------------------*/ -class dsmcCloud +template +class MixedDiffuseSpecular : - public DsmcCloud + public WallInteractionModel { - // Private member functions - - //- Disallow default bitwise copy construct - dsmcCloud(const dsmcCloud&); + // Private data - //- Disallow default bitwise assignment - void operator=(const dsmcCloud&); + //- Fraction of wall interactions that are diffuse + scalar diffuseFraction_; public: //- Runtime type information - TypeName("dsmcCloud"); + TypeName("MixedDiffuseSpecular"); // Constructors - //- Construct from components - dsmcCloud + //- Construct from dictionary + MixedDiffuseSpecular ( - const word& cloudName, - const volScalarField& T, - const volVectorField& U + const dictionary& dict, + CloudType& cloud ); - //- Construct from name and mesh, used to initialise. - dsmcCloud - ( - const word& cloudName, - const fvMesh& mesh - ); - //- Destructor - ~dsmcCloud(); + // Destructor + virtual ~MixedDiffuseSpecular(); - // Member functions + // Member Functions - //- Write fields - virtual void writeFields() const; + //- 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 // ************************************************************************* // diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C index 8957dee1..cbc00df7 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; @@ -178,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; @@ -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 -- 2.11.4.GIT