diff --git a/src/lagrangian/dsmc/Make/files b/src/lagrangian/dsmc/Make/files index 2359f74dec..9f6d33d209 100644 --- a/src/lagrangian/dsmc/Make/files +++ b/src/lagrangian/dsmc/Make/files @@ -11,5 +11,6 @@ clouds/derived/dsmcCloud/dsmcCloud.C parcels/derived/dsmcParcel/defineDsmcParcel.C parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C +parcels/derived/dsmcParcel/makeDsmcParcelInflowBoundaryModels.C LIB = $(FOAM_LIBBIN)/libdsmc diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C index 1867670f70..b6d20bd630 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C @@ -27,6 +27,7 @@ License #include "DsmcCloud.H" #include "BinaryCollisionModel.H" #include "WallInteractionModel.H" +#include "InflowBoundaryModel.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -234,8 +235,8 @@ void Foam::DsmcCloud::collisions() scalar sigmaTcRMax = sigmaTcRMax_[celli]; scalar selectedPairs = collisionSelectionRemainder_[celli] - + 0.5*nC*(nC-1)*nParticle_*sigmaTcRMax*deltaT - /mesh_.cellVolumes()[celli]; + + 0.5*nC*(nC-1)*nParticle_*sigmaTcRMax*deltaT + /mesh_.cellVolumes()[celli]; label nCandidates(selectedPairs); @@ -403,6 +404,14 @@ Foam::DsmcCloud::DsmcCloud particleProperties_, *this ) + ), + inflowBoundaryModel_ + ( + InflowBoundaryModel >::New + ( + particleProperties_, + *this + ) ) { buildConstProps(); @@ -490,7 +499,8 @@ Foam::DsmcCloud::DsmcCloud ) ), binaryCollisionModel_(), - wallInteractionModel_() + wallInteractionModel_(), + inflowBoundaryModel_() { clear(); @@ -526,13 +536,14 @@ void Foam::DsmcCloud::evolve() { typename ParcelType::trackData td(*this); - //this->injection().inject(td); - if (debug) { this->dumpParticlePositions(); } + // Insert new particles from the inflow boundary + this->inflowBoundary().inflow(); + // Move the particles ballistically with their current velocities Cloud::move(td); diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H index 8ea6f31fe1..5d9edc9c06 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H @@ -59,6 +59,9 @@ class BinaryCollisionModel; template class WallInteractionModel; +template +class InflowBoundaryModel; + /*---------------------------------------------------------------------------*\ Class DsmcCloud Declaration \*---------------------------------------------------------------------------*/ @@ -127,6 +130,10 @@ class DsmcCloud autoPtr > > wallInteractionModel_; + //- Inflow boundary model + autoPtr > > + inflowBoundaryModel_; + // Private Member Functions @@ -205,7 +212,7 @@ public: //- Return the number of real particles represented by one // parcel - inline label nParticle() const; + inline scalar nParticle() const; //- Return the cell occupancy addressing inline const List >& @@ -292,6 +299,14 @@ public: inline WallInteractionModel >& wallInteraction(); + //- Return reference to wall interaction model + inline const InflowBoundaryModel >& + inflowBoundary() const; + + //- Return non-const reference to wall interaction model + inline InflowBoundaryModel >& + inflowBoundary(); + // Check diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H index 057570030e..39c25f5227 100644 --- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H +++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H @@ -57,7 +57,7 @@ Foam::DsmcCloud::typeIdList() const template -inline Foam::label Foam::DsmcCloud::nParticle() const +inline Foam::scalar Foam::DsmcCloud::nParticle() const { return nParticle_; } @@ -166,6 +166,22 @@ Foam::DsmcCloud::wallInteraction() } +template +inline const Foam::InflowBoundaryModel >& +Foam::DsmcCloud::inflowBoundary() const +{ + return inflowBoundaryModel_; +} + + +template +inline Foam::InflowBoundaryModel >& +Foam::DsmcCloud::inflowBoundary() +{ + return inflowBoundaryModel_(); +} + + template inline Foam::scalar Foam::DsmcCloud::massInSystem() const { diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelInflowBoundaryModels.C b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelInflowBoundaryModels.C new file mode 100644 index 0000000000..c3ee121ff0 --- /dev/null +++ b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelInflowBoundaryModels.C @@ -0,0 +1,52 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-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 "dsmcParcel.H" +#include "DsmcCloud.H" +#include "FreeStream.H" +#include "NoInflow.H" + +namespace Foam +{ + makeInflowBoundaryModel(DsmcCloud); + + // Add instances of inflow boundary model to the table + makeInflowBoundaryModelType + ( + FreeStream, + DsmcCloud, + dsmcParcel + ); + makeInflowBoundaryModelType + ( + NoInflow, + DsmcCloud, + dsmcParcel + ); +}; + + +// ************************************************************************* // diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C index 5f00f4ca93..1cf6fd90dc 100644 --- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C +++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.C @@ -109,7 +109,8 @@ Foam::LarsenBorgnakkeVariableHardSphere::LarsenBorgnakkeVariableHardS // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template -Foam::LarsenBorgnakkeVariableHardSphere::~LarsenBorgnakkeVariableHardSphere() +Foam::LarsenBorgnakkeVariableHardSphere:: +~LarsenBorgnakkeVariableHardSphere() {} diff --git a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.H b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.H index 1e0f8bd975..c95d70f24c 100644 --- a/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.H +++ b/src/lagrangian/dsmc/submodels/BinaryCollisionModel/LarsenBorgnakkeVariableHardSphere/LarsenBorgnakkeVariableHardSphere.H @@ -50,7 +50,7 @@ class LarsenBorgnakkeVariableHardSphere { // Private data - //- Temperature + //- Relaxation collision number const scalar relaxationCollisionNumber_; diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C new file mode 100644 index 0000000000..ed55544077 --- /dev/null +++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C @@ -0,0 +1,218 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-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 "FreeStream.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::FreeStream::FreeStream +( + const dictionary& dict, + CloudType& cloud +) +: + InflowBoundaryModel(dict, cloud, typeName), + patchIndex_(), + temperature_(readScalar(this->coeffDict().lookup("temperature"))), + velocity_(this->coeffDict().lookup("velocity")), + moleculeTypeIds_(), + numberDensities_(), + particleFluxAccumulators_() +{ + word patchName = this->coeffDict().lookup("patch"); + + patchIndex_ = cloud.mesh().boundaryMesh().findPatchID(patchName); + + const polyPatch& patch = cloud.mesh().boundaryMesh()[patchIndex_]; + + if (patchIndex_ == -1) + { + FatalErrorIn("Foam::DsmcCloud::initialise") + << "patch " << patchName << " not found." << nl + << abort(FatalError); + } + + const dictionary& numberDensitiesDict + ( + this->coeffDict().subDict("numberDensities") + ); + + List molecules(numberDensitiesDict.toc()); + + numberDensities_.setSize(molecules.size()); + + moleculeTypeIds_.setSize(molecules.size()); + + forAll(molecules, i) + { + numberDensities_[i] = readScalar + ( + numberDensitiesDict.lookup(molecules[i]) + ); + + moleculeTypeIds_[i] = findIndex(cloud.typeIdList(), molecules[i]); + + if (moleculeTypeIds_[i] == -1) + { + FatalErrorIn("Foam::DsmcCloud::initialise") + << "typeId " << molecules[i] << "not defined in cloud." << nl + << abort(FatalError); + } + } + + numberDensities_ /= cloud.nParticle(); + + particleFluxAccumulators_.setSize + ( + molecules.size(), + Field(patch.size(), 0) + ); +} + + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::FreeStream::~FreeStream() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::FreeStream::inflow() +{ + CloudType& cloud(this->owner()); + + const polyMesh& mesh(cloud.mesh()); + + const scalar deltaT = mesh.time().deltaT().value(); + + Random& rndGen(cloud.rndGen()); + + const polyPatch& patch = mesh.boundaryMesh()[patchIndex_]; + + label particlesInserted = 0; + + // Add mass to the accumulators. negative face area dotted with the + // velocity to point flux into the domain. + + forAll(particleFluxAccumulators_, i) + { + particleFluxAccumulators_[i] += + -patch.faceAreas() & (velocity_*numberDensities_[i]*deltaT); + + forAll(particleFluxAccumulators_[i], f) + { + scalar& faceAccumulator = particleFluxAccumulators_[i][f]; + + // Number of particles to insert + label nI = max(label(faceAccumulator), 0); + + faceAccumulator -= nI; + + label typeId = moleculeTypeIds_[i]; + + scalar mass = cloud.constProps(typeId).mass(); + + labelList faceVertices = patch[f]; + + label globalFaceIndex = f + patch.start(); + + label cell = mesh.faceOwner()[globalFaceIndex]; + + const vector& fC = patch.faceCentres()[f]; + + for (label n = 0; n < nI; n++) + { + // Temporarily insert particles half way between the face and + // cell centres + vector p = 0.5*(fC + mesh.cellCentres()[cell]); + + // Normal unit vector *negative* so normal is pointing into the + // domain + vector nw = patch.faceAreas()[f]; + nw /= -mag(nw); + + // Wall tangential unit vector. Use the direction between the + // face centre and the first vertex in the list + vector tw1 = fC - (mesh.points()[faceVertices[0]]); + tw1 /= mag(tw1); + + // Other tangential unit vector. Rescaling in case face is not + // flat and nw and tw1 aren't perfectly orthogonal + vector tw2 = nw ^ tw1; + tw2 /= mag(tw2); + + scalar C = sqrt(CloudType::kb*temperature_/mass); + + vector U = + C + *( + rndGen.GaussNormal()*tw1 + + rndGen.GaussNormal()*tw2 + - sqrt(-2.0*log(max(1 - rndGen.scalar01(),VSMALL)))*nw + ); + + U += velocity_; + + scalar Ei = + 0.5*cloud.constProps(typeId).internalDegreesOfFreedom() + *CloudType::kb*temperature_; + + cloud.addNewParcel + ( + p, + U, + Ei, + cell, + typeId + ); + + particlesInserted++; + } + } + } + + reduce(particlesInserted, sumOp