FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

FreeStream.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "FreeStream.H"
00027 
00028 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00029 
00030 template <class CloudType>
00031 Foam::FreeStream<CloudType>::FreeStream
00032 (
00033     const dictionary& dict,
00034     CloudType& cloud
00035 )
00036 :
00037     InflowBoundaryModel<CloudType>(dict, cloud, typeName),
00038     patches_(),
00039     moleculeTypeIds_(),
00040     numberDensities_(),
00041     particleFluxAccumulators_()
00042 {
00043     // Identify which patches to use
00044 
00045     DynamicList<label> patches;
00046 
00047     forAll(cloud.mesh().boundaryMesh(), p)
00048     {
00049         const polyPatch& patch = cloud.mesh().boundaryMesh()[p];
00050 
00051         if (isType<polyPatch>(patch))
00052         {
00053             patches.append(p);
00054         }
00055     }
00056 
00057     patches_.transfer(patches);
00058 
00059     const dictionary& numberDensitiesDict
00060     (
00061         this->coeffDict().subDict("numberDensities")
00062     );
00063 
00064     List<word> molecules(numberDensitiesDict.toc());
00065 
00066     // Initialise the particleFluxAccumulators_
00067     particleFluxAccumulators_.setSize(patches_.size());
00068 
00069     forAll(patches_, p)
00070     {
00071         const polyPatch& patch = cloud.mesh().boundaryMesh()[patches_[p]];
00072 
00073         particleFluxAccumulators_[p] = List<Field<scalar> >
00074         (
00075             molecules.size(),
00076             Field<scalar>(patch.size(), 0.0)
00077         );
00078     }
00079 
00080     moleculeTypeIds_.setSize(molecules.size());
00081 
00082     numberDensities_.setSize(molecules.size());
00083 
00084     forAll(molecules, i)
00085     {
00086         numberDensities_[i] = readScalar
00087         (
00088             numberDensitiesDict.lookup(molecules[i])
00089         );
00090 
00091         moleculeTypeIds_[i] = findIndex(cloud.typeIdList(), molecules[i]);
00092 
00093         if (moleculeTypeIds_[i] == -1)
00094         {
00095             FatalErrorIn
00096             (
00097                 "Foam::FreeStream<CloudType>::FreeStream"
00098                 "("
00099                     "const dictionary&, "
00100                     "CloudType&"
00101                 ")"
00102             )   << "typeId " << molecules[i] << "not defined in cloud." << nl
00103                 << abort(FatalError);
00104         }
00105     }
00106 
00107     numberDensities_ /= cloud.nParticle();
00108 }
00109 
00110 
00111 
00112 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00113 
00114 template <class CloudType>
00115 Foam::FreeStream<CloudType>::~FreeStream()
00116 {}
00117 
00118 
00119 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00120 
00121 template <class CloudType>
00122 void Foam::FreeStream<CloudType>::inflow()
00123 {
00124     CloudType& cloud(this->owner());
00125 
00126     const polyMesh& mesh(cloud.mesh());
00127 
00128     const scalar deltaT = mesh.time().deltaTValue();
00129 
00130     Random& rndGen(cloud.rndGen());
00131 
00132     scalar sqrtPi = sqrt(mathematicalConstant::pi);
00133 
00134     label particlesInserted = 0;
00135 
00136     const volScalarField::GeometricBoundaryField& boundaryT
00137     (
00138         cloud.boundaryT().boundaryField()
00139     );
00140 
00141     const volVectorField::GeometricBoundaryField& boundaryU
00142     (
00143         cloud.boundaryU().boundaryField()
00144     );
00145 
00146     forAll(patches_, p)
00147     {
00148         label patchI = patches_[p];
00149 
00150         const polyPatch& patch = mesh.boundaryMesh()[patchI];
00151 
00152         // Add mass to the accumulators.  negative face area dotted with the
00153         // velocity to point flux into the domain.
00154 
00155         // Take a reference to the particleFluxAccumulator for this patch
00156         List<Field<scalar> >& pFA = particleFluxAccumulators_[p];
00157 
00158         forAll(pFA, i)
00159         {
00160             label typeId = moleculeTypeIds_[i];
00161 
00162             scalar mass = cloud.constProps(typeId).mass();
00163 
00164             if (min(boundaryT[patchI]) < SMALL)
00165             {
00166                 FatalErrorIn ("Foam::FreeStream<CloudType>::inflow()")
00167                     << "Zero boundary temperature detected, "
00168                     << "check boundaryT condition." << nl
00169                     << nl << abort(FatalError);
00170             }
00171 
00172             scalarField mostProbableSpeed
00173             (
00174                 cloud.maxwellianMostProbableSpeed
00175                 (
00176                     boundaryT[patchI],
00177                     mass
00178                 )
00179             );
00180 
00181             // Dotting boundary velocity with the face unit normal (which points
00182             // out of the domain, so it must be negated), dividing by the most
00183             // probable speed to form molecularSpeedRatio * cosTheta
00184 
00185             scalarField sCosTheta =
00186                 boundaryU[patchI]
00187               & -patch.faceAreas()/mag(patch.faceAreas())
00188                /mostProbableSpeed;
00189 
00190             // From Bird eqn 4.22
00191 
00192             pFA[i] +=
00193                 mag(patch.faceAreas())*numberDensities_[i]*deltaT
00194                *mostProbableSpeed
00195                *(
00196                    exp(-sqr(sCosTheta)) + sqrtPi*sCosTheta*(1 + erf(sCosTheta))
00197                 )
00198                /(2.0*sqrtPi);
00199         }
00200 
00201         forAll(patch, f)
00202         {
00203             // Loop over all faces as the outer loop to avoid calculating
00204             // geometrical properties multiple times.
00205 
00206             labelList faceVertices = patch[f];
00207 
00208             label nVertices = faceVertices.size();
00209 
00210             label globalFaceIndex = f + patch.start();
00211 
00212             label cell = mesh.faceOwner()[globalFaceIndex];
00213 
00214             const vector& fC = patch.faceCentres()[f];
00215 
00216             scalar fA = mag(patch.faceAreas()[f]);
00217 
00218             // Cumulative triangle area fractions
00219             List<scalar> cTriAFracs(nVertices);
00220             cTriAFracs[0] = 0.0;
00221 
00222             for (label v = 0; v < nVertices - 1; v++)
00223             {
00224                 const point& vA = mesh.points()[faceVertices[v]];
00225 
00226                 const point& vB = mesh.points()[faceVertices[(v + 1)]];
00227 
00228                 cTriAFracs[v] =
00229                     0.5*mag((vA - fC)^(vB - fC))/fA
00230                   + cTriAFracs[max((v - 1), 0)];
00231             }
00232 
00233             // Force the last area fraction value to 1.0 to avoid any
00234             // rounding/non-flat face errors giving a value < 1.0
00235             cTriAFracs[nVertices - 1] = 1.0;
00236 
00237             // Normal unit vector *negative* so normal is pointing into the
00238             // domain
00239             vector n = patch.faceAreas()[f];
00240             n /= -mag(n);
00241 
00242             // Wall tangential unit vector. Use the direction between the
00243             // face centre and the first vertex in the list
00244             vector t1 = fC - (mesh.points()[faceVertices[0]]);
00245             t1 /= mag(t1);
00246 
00247             // Other tangential unit vector.  Rescaling in case face is not
00248             // flat and n and t1 aren't perfectly orthogonal
00249             vector t2 = n^t1;
00250             t2 /= mag(t2);
00251 
00252             scalar faceTemperature = boundaryT[patchI][f];
00253 
00254             const vector& faceVelocity = boundaryU[patchI][f];
00255 
00256             forAll(pFA, i)
00257             {
00258                 scalar& faceAccumulator = pFA[i][f];
00259 
00260                 // Number of whole particles to insert
00261                 label nI = max(label(faceAccumulator), 0);
00262 
00263                 // Add another particle with a probability proportional to the
00264                 // remainder of taking the integer part of faceAccumulator
00265                 if ((faceAccumulator - nI) > rndGen.scalar01())
00266                 {
00267                     nI++;
00268                 }
00269 
00270                 faceAccumulator -= nI;
00271 
00272                 label typeId = moleculeTypeIds_[i];
00273 
00274                 scalar mass = cloud.constProps(typeId).mass();
00275 
00276                 for (label i = 0; i < nI; i++)
00277                 {
00278                     // Choose a triangle to insert on, based on their relative
00279                     // area
00280 
00281                     scalar triSelection = rndGen.scalar01();
00282 
00283                     // Selected triangle
00284                     label sTri = -1;
00285 
00286                     forAll(cTriAFracs, tri)
00287                     {
00288                         sTri = tri;
00289 
00290                         if (cTriAFracs[tri] >= triSelection)
00291                         {
00292                             break;
00293                         }
00294                     }
00295 
00296                     // Randomly distribute the points on the triangle, using the
00297                     // method from:
00298                     // Generating Random Points in Triangles
00299                     // by Greg Turk
00300                     // from "Graphics Gems", Academic Press, 1990
00301                     // http://tog.acm.org/GraphicsGems/gems/TriPoints.c
00302 
00303                     const point& A = fC;
00304                     const point& B = mesh.points()[faceVertices[sTri]];
00305                     const point& C =
00306                     mesh.points()[faceVertices[(sTri + 1) % nVertices]];
00307 
00308                     scalar s = rndGen.scalar01();
00309                     scalar t = sqrt(rndGen.scalar01());
00310 
00311                     point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
00312 
00313                     // Velocity generation
00314 
00315                     scalar mostProbableSpeed
00316                     (
00317                         cloud.maxwellianMostProbableSpeed
00318                         (
00319                             faceTemperature,
00320                             mass
00321                         )
00322                     );
00323 
00324                     scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
00325 
00326                     // Coefficients required for Bird eqn 12.5
00327                     scalar uNormProbCoeffA =
00328                         sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
00329 
00330                     scalar uNormProbCoeffB =
00331                         0.5*
00332                         (
00333                             1.0
00334                           + sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
00335                         );
00336 
00337                     // Equivalent to the QA value in Bird's DSMC3.FOR
00338                     scalar randomScaling = 3.0;
00339 
00340                     if (sCosTheta < -3)
00341                     {
00342                         randomScaling = mag(sCosTheta) + 1;
00343                     }
00344 
00345                     scalar P = -1;
00346 
00347                     // Normalised candidates for the normal direction velocity
00348                     // component
00349                     scalar uNormal;
00350                     scalar uNormalThermal;
00351 
00352                     // Select a velocity using Bird eqn 12.5
00353                     do
00354                     {
00355                         uNormalThermal =
00356                             randomScaling*(2.0*rndGen.scalar01() - 1);
00357 
00358                         uNormal = uNormalThermal + sCosTheta;
00359 
00360                         if (uNormal < 0.0)
00361                         {
00362                             P = -1;
00363                         }
00364                         else
00365                         {
00366                             P = 2.0*uNormal/uNormProbCoeffA
00367                                *exp(uNormProbCoeffB - sqr(uNormalThermal));
00368                         }
00369 
00370                     } while (P < rndGen.scalar01());
00371 
00372                     vector U =
00373                         sqrt(CloudType::kb*faceTemperature/mass)
00374                        *(
00375                             rndGen.GaussNormal()*t1
00376                           + rndGen.GaussNormal()*t2
00377                         )
00378                       + (t1 & faceVelocity)*t1
00379                       + (t2 & faceVelocity)*t2
00380                       + mostProbableSpeed*uNormal*n;
00381 
00382                     scalar Ei = cloud.equipartitionInternalEnergy
00383                     (
00384                         faceTemperature,
00385                         cloud.constProps(typeId).internalDegreesOfFreedom()
00386                     );
00387 
00388                     cloud.addNewParcel
00389                     (
00390                         p,
00391                         U,
00392                         Ei,
00393                         cell,
00394                         typeId
00395                     );
00396 
00397                     particlesInserted++;
00398                 }
00399             }
00400         }
00401     }
00402 
00403     reduce(particlesInserted, sumOp<label>());
00404 
00405     Info<< "    Particles inserted              = "
00406         << particlesInserted << endl;
00407 
00408 }
00409 
00410 
00411 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines