/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ 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 3 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, see . \*---------------------------------------------------------------------------*/ #include "DSMCCloud.H" #include "BinaryCollisionModel.H" #include "WallInteractionModel.H" #include "InflowBoundaryModel.H" #include "constants.H" #include "zeroGradientFvPatchFields.H" #include "polyMeshTetDecomposition.H" using namespace Foam::constant; // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template void Foam::DSMCCloud::buildConstProps() { Info<< nl << "Constructing constant properties for" << endl; constProps_.setSize(typeIdList_.size()); dictionary moleculeProperties ( particleProperties_.subDict("moleculeProperties") ); forAll(typeIdList_, i) { const word& id(typeIdList_[i]); Info<< " " << id << endl; const dictionary& molDict(moleculeProperties.subDict(id)); constProps_[i] = typename ParcelType::constantProperties(molDict); } } template void Foam::DSMCCloud::buildCellOccupancy() { forAll(cellOccupancy_, cO) { cellOccupancy_[cO].clear(); } forAllIter(typename DSMCCloud, *this, iter) { cellOccupancy_[iter().cell()].append(&iter()); } } template void Foam::DSMCCloud::initialise ( const IOdictionary& dsmcInitialiseDict ) { Info<< nl << "Initialising particles" << endl; const scalar temperature ( readScalar(dsmcInitialiseDict.lookup("temperature")) ); const vector velocity(dsmcInitialiseDict.lookup("velocity")); const dictionary& numberDensitiesDict ( dsmcInitialiseDict.subDict("numberDensities") ); List molecules(numberDensitiesDict.toc()); Field numberDensities(molecules.size()); forAll(molecules, i) { numberDensities[i] = readScalar ( numberDensitiesDict.lookup(molecules[i]) ); } numberDensities /= nParticle_; forAll(mesh_.cells(), celli) { List cellTets = polyMeshTetDecomposition::cellTetIndices ( mesh_, celli ); forAll(cellTets, tetI) { const tetIndices& cellTetIs = cellTets[tetI]; tetPointRef tet = cellTetIs.tet(mesh_); scalar tetVolume = tet.mag(); forAll(molecules, i) { const word& moleculeName(molecules[i]); label typeId = typeIdList_.find(moleculeName); if (typeId == -1) { FatalErrorInFunction << "typeId " << moleculeName << "not defined." << nl << abort(FatalError); } const typename ParcelType::constantProperties& cP = constProps(typeId); scalar numberDensity = numberDensities[i]; // Calculate the number of particles required scalar particlesRequired = numberDensity*tetVolume; // Only integer numbers of particles can be inserted label nParticlesToInsert = label(particlesRequired); // Add another particle with a probability proportional to the // remainder of taking the integer part of particlesRequired if ( (particlesRequired - nParticlesToInsert) > rndGen_.sample01() ) { nParticlesToInsert++; } for (label pI = 0; pI < nParticlesToInsert; pI++) { point p = tet.randomPoint(rndGen_); vector U = equipartitionLinearVelocity ( temperature, cP.mass() ); scalar Ei = equipartitionInternalEnergy ( temperature, cP.internalDegreesOfFreedom() ); U += velocity; addNewParcel(p, celli, U, Ei, typeId); } } } } // Initialise the sigmaTcRMax_ field to the product of the cross section of // the most abundant species and the most probable thermal speed (Bird, // p222-223) label mostAbundantType(findMax(numberDensities)); const typename ParcelType::constantProperties& cP = constProps ( mostAbundantType ); sigmaTcRMax_.primitiveFieldRef() = cP.sigmaT()*maxwellianMostProbableSpeed ( temperature, cP.mass() ); sigmaTcRMax_.correctBoundaryConditions(); } template void Foam::DSMCCloud::collisions() { if (!binaryCollision().active()) { return; } // Temporary storage for subCells List> subCells(8); scalar deltaT = mesh().time().deltaTValue(); label collisionCandidates = 0; label collisions = 0; forAll(cellOccupancy_, celli) { const DynamicList& cellParcels(cellOccupancy_[celli]); label nC(cellParcels.size()); if (nC > 1) { // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // Assign particles to one of 8 Cartesian subCells // Clear temporary lists forAll(subCells, i) { subCells[i].clear(); } // Inverse addressing specifying which subCell a parcel is in List