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 84b499451f..f5a57f15e8 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 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/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/etc/controlDict b/etc/controlDict index 97af349e7a..ef72407eaf 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/dsmc/Make/files b/src/lagrangian/dsmc/Make/files index 9f6d33d209..983de123ab 100644 --- a/src/lagrangian/dsmc/Make/files +++ b/src/lagrangian/dsmc/Make/files @@ -4,9 +4,6 @@ parcels/derived/dsmcParcel/dsmcParcel.C /* Cloud base classes */ clouds/baseClasses/DsmcBaseCloud/DsmcBaseCloud.C -/* Clouds */ -clouds/derived/dsmcCloud/dsmcCloud.C - /* submodels */ parcels/derived/dsmcParcel/defineDsmcParcel.C parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C diff --git a/src/lagrangian/dsmc/Make/options b/src/lagrangian/dsmc/Make/options index 15874b7b55..55865dfabc 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 5ac29101e5..f7aaaa724d 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C @@ -297,7 +297,7 @@ void Foam::DsmcCloud::collisions() // Temporary storage for subCells List > subCells(8); - scalar deltaT = cachedDeltaT(); + scalar deltaT = mesh().time().deltaTValue(); label collisionCandidates = 0; @@ -473,21 +473,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); - forAll(qBF, p) + 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 + ); +} + + +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(); + + forAllConstIter(typename DsmcCloud, *this, iter) { - qBF[p] = 0.0; + 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(); } - volVectorField::GeometricBoundaryField& fDBF(fD_.boundaryField()); + rhoN *= nParticle_/mesh().cellVolumes(); + rhoN_.correctBoundaryConditions(); - forAll(fDBF, p) - { - fDBF[p] = vector::zero; - } + 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(); } @@ -523,15 +594,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 @@ -563,37 +633,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 @@ -641,7 +816,8 @@ template Foam::DsmcCloud::DsmcCloud ( const word& cloudName, - const fvMesh& mesh + const fvMesh& mesh, + const IOdictionary& dsmcInitialiseDict ) : Cloud(mesh, cloudName, false), @@ -707,15 +883,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, @@ -725,13 +997,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, @@ -754,18 +1026,6 @@ Foam::DsmcCloud::DsmcCloud buildConstProps(); - IOdictionary dsmcInitialiseDict - ( - IOobject - ( - "dsmcInitialiseDict", - mesh_.time().system(), - mesh_, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ); - initialise(dsmcInitialiseDict); } @@ -782,13 +1042,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) { @@ -803,6 +1060,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 c3765ba0df..5b2f382eb6 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&); @@ -173,20 +193,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 ); @@ -242,35 +263,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 force density at surface field - inline volVectorField& fD(); + //- 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(); // 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 @@ -385,29 +443,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 diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H index 9b65d77659..beaacb2276 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H @@ -126,58 +126,82 @@ inline Foam::Random& Foam::DsmcCloud::rndGen() template -inline void Foam::DsmcCloud::storeDeltaT() +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::qBF() { - cachedDeltaT_ = mesh().time().deltaTValue(); + 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 const Foam::volScalarField& Foam::DsmcCloud::T() const +inline Foam::volScalarField::GeometricBoundaryField& +Foam::DsmcCloud::iDofBF() { - return T_; + return iDof_.boundaryField(); } template -inline const Foam::volVectorField& Foam::DsmcCloud::U() const +inline Foam::volVectorField::GeometricBoundaryField& +Foam::DsmcCloud::momentumBF() { - return U_; + return momentum_.boundaryField(); +} + + +template +inline const Foam::volScalarField& +Foam::DsmcCloud::boundaryT() const +{ + return boundaryT_; +} + + +template +inline const Foam::volVectorField& +Foam::DsmcCloud::boundaryU() const +{ + return boundaryU_; } @@ -381,265 +405,70 @@ Foam::DsmcCloud::maxwellianMostProbableSpeed template -inline const Foam::tmp +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::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::tmp -Foam::DsmcCloud::rhoM() const +inline const Foam::volScalarField& Foam::DsmcCloud::rhoM() 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 rhoM_; } 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) - ) - ); - - scalarField& dsmcRhoN = tdsmcRhoN().internalField(); - forAllConstIter(typename DsmcCloud, *this, iter) - { - const ParcelType& p = iter(); - const label cellI = p.cell(); - - dsmcRhoN[cellI]++; - } - - return tdsmcRhoN; + return dsmcRhoN_; } template -inline const Foam::tmp -Foam::DsmcCloud::momentum() const -{ - tmp tmomentum - ( - new volVectorField - ( - IOobject - ( - this->name() + "momentum", - this->db().time().timeName(), - this->db(), - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh_, - dimensionedVector - ( - "zero", - dimensionSet(1, -2, -1, 0, 0), - vector::zero - ) - ) - ); - - vectorField& momentum = tmomentum().internalField(); - forAllConstIter(typename DsmcCloud, *this, iter) - { - const ParcelType& p = iter(); - const label cellI = p.cell(); - - momentum[cellI] += constProps(p.typeId()).mass()*p.U(); - } - - momentum *= nParticle_/mesh().cellVolumes(); - - return tmomentum; -} - - -template -inline const Foam::tmp +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) - ) - ); + return iDof_; +} - 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::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 deleted file mode 100644 index 7f94f3e739..0000000000 --- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C +++ /dev/null @@ -1,66 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 "dsmcCloud.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ - defineTypeNameAndDebug(dsmcCloud, 0); -}; - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::dsmcCloud::dsmcCloud -( - const word& cloudName, - const volScalarField& T, - const volVectorField& U -) -: - DsmcCloud(cloudName, T, U) -{} - - -Foam::dsmcCloud::dsmcCloud -( - const word& cloudName, - const fvMesh& mesh -) -: - DsmcCloud(cloudName, mesh) -{} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::dsmcCloud::~dsmcCloud() -{} - - -// ************************************************************************* // diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H index 99343c2b3e..c0c17336e6 100644 --- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H +++ b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H @@ -43,54 +43,8 @@ SourceFiles namespace Foam { - -/*---------------------------------------------------------------------------*\ - Class dsmcCloud Declaration -\*---------------------------------------------------------------------------*/ - -class dsmcCloud -: - public DsmcCloud -{ - // Private member functions - - //- Disallow default bitwise copy construct - dsmcCloud(const dsmcCloud&); - - //- Disallow default bitwise assignment - void operator=(const dsmcCloud&); - - -public: - - //- Runtime type information - TypeName("dsmcCloud"); - - - // Constructors - - //- Construct from components - dsmcCloud - ( - const word& cloudName, - const volScalarField& T, - const volVectorField& U - ); - - //- Construct from name and mesh, used to initialise. - dsmcCloud - ( - const word& cloudName, - const fvMesh& mesh - ); - - //- Destructor - ~dsmcCloud(); -}; - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam + typedef DsmcCloud dsmcCloud; +} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C index 521559b72c..3d752aa98d 100644 --- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C +++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C @@ -25,6 +25,7 @@ License \*---------------------------------------------------------------------------*/ #include "DsmcParcel.H" +#include "meshTools.H" // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // @@ -43,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; @@ -113,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_; @@ -132,27 +179,40 @@ void Foam::DsmcParcel::hitWallPatch typeId_ ); + U_dot_nw = U_ & nw; + + Ut = U_ - U_dot_nw*nw; + + 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; + // 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; } @@ -212,4 +272,3 @@ void Foam::DsmcParcel::transformProperties #include "DsmcParcelIO.C" // ************************************************************************* // - 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/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C index 03fc89d032..1df22c68d5 100644 --- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C +++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C @@ -129,7 +129,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()); @@ -139,12 +139,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) @@ -168,7 +168,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 569e9b21e6..38890400c8 100644 --- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C +++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C @@ -71,10 +71,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()); @@ -93,9 +93,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 @@ -104,7 +104,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(); @@ -118,7 +118,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 new file mode 100644 index 0000000000..9f2f5fdc3b --- /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 U_dot_nw = U & nw; + + CloudType& cloud(this->owner()); + + Random& rndGen(cloud.rndGen()); + + if (diffuseFraction_ > rndGen.scalar01()) + { + // Diffuse reflection + + // Wall tangential velocity (flow direction) + vector Ut = U - U_dot_nw*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()) + ); + + U_dot_nw = U & nw; + + Ut = U - U_dot_nw*nw; + } + + // Wall tangential unit vector + vector tw1 = Ut/mag(Ut); + + // Other tangential unit vector + vector tw2 = nw^tw1; + + scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace]; + + scalar mass = cloud.constProps(typeId).mass(); + + scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom(); + + U = + sqrt(physicoChemical::k.value()*T/mass) + *( + rndGen.GaussNormal()*tw1 + + rndGen.GaussNormal()*tw2 + - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw + ); + + U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace]; + + Ei = cloud.equipartitionInternalEnergy(T, iDof); + } + else + { + // Specular reflection + + if (U_dot_nw > 0.0) + { + U -= 2.0*U_dot_nw*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 + +// ************************************************************************* // 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; } } 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