zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge commit 'OpenCFD/master' into olesenm

Browse Source
This commit is contained in:
Mark Olesen
2009-11-23 12:03:10 +01:00
parent 534b7dca70 e06171fea4
commit 87300f124e
49 changed files with 1070 additions and 962 deletions
Download Patch File Download Diff File Expand all files Collapse all files

162
applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H
View File

@ -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);

42
applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
View File

@ -41,53 +41,21 @@ int main(int argc, char *argv[])
#include "setRootCase.H" #include "setRootCase.H"
#include "createTime.H" #include "createTime.H"
#include "createMesh.H" #include "createMesh.H"
#include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< nl << "Constructing dsmcCloud " << endl;
dsmcCloud dsmc("dsmc", mesh);
Info<< "\nStarting time loop\n" << endl; Info<< "\nStarting time loop\n" << endl;
while (runTime.run()) while (runTime.loop())
{ {
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl; Info<< "Time = " << runTime.timeName() << nl << endl;
// Carry out dsmcCloud timestep
dsmc.evolve(); 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(); dsmc.info();
runTime.write(); runTime.write();

133
applications/test/DimensionedField/DimensionedFieldTest.C
View File

@ -1,133 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "argList.H"
#include "Time.H"
#include "DimensionedFields.H"
#include "DimensionedSphericalTensorField.H"
#include "vector.H"
#include "tensor.H"
#include "GeoMesh.H"
using namespace Foam;
namespace Foam
{
class vMesh
{
public:
vMesh()
{}
label size() const
{
return 10;
}
};
};
template<>
const word Foam::DimensionedField<scalar, GeoMesh<vMesh> >::typeName
(
"dimenionedScalarField"
);
template<>
const word Foam::DimensionedField<vector, GeoMesh<vMesh> >::typeName
(
"dimenionedVectorField"
);
template<>
const word Foam::DimensionedField<tensor, GeoMesh<vMesh> >::typeName
(
"dimenionedTensorField"
);
template<>
const word Foam::DimensionedField<sphericalTensor, GeoMesh<vMesh> >::typeName
(
"dimenionedSphericalTensorField"
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
vMesh vm;
DimensionedField<scalar, GeoMesh<vMesh> > dsf
(
IOobject
(
"dsf",
runTime.timeName(),
runTime,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
vm
);
Info<< dsf << endl;
dsf += dsf;
dsf -= dimensionedScalar("5", dsf.dimensions(), 5.0);
Info<< dsf << endl;
Info<< sqr(dsf + dsf) - sqr(dsf + dsf) << endl;
DimensionedField<vector, GeoMesh<vMesh> > dvf
(
IOobject
(
"dvf",
runTime.timeName(),
runTime,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
vm
);
Info<< (dvf ^ (dvf ^ dvf)) << endl;
Info << "End\n" << endl;
return 0;
}
// ************************************************************************* //

3
applications/test/DimensionedField/Make/files
View File

@ -1,3 +0,0 @@
DimensionedFieldTest.C
EXE = $(FOAM_USER_APPBIN)/DimensionedFieldTest

2
applications/test/DimensionedField/Make/options
View File

@ -1,2 +0,0 @@
/* EXE_INC = -I$(LIB_SRC)/cfdTools/include */
/* EXE_LIBS = -lfiniteVolume */

3
applications/test/HashTable3/hashTableTest3.C
View File

@ -33,7 +33,6 @@ Description
#include "HashPtrTable.H" #include "HashPtrTable.H"
#include "Map.H" #include "Map.H"
#include "StaticHashTable.H" #include "StaticHashTable.H"
#include "HashTbl.H"
#include "cpuTime.H" #include "cpuTime.H"
using namespace Foam; using namespace Foam;
@ -53,7 +52,7 @@ int main(int argc, char *argv[])
// Map<label> map(2 * nSize); // Map<label> map(2 * nSize);
// HashTable<label, label, Hash<label> > map(2 * nSize); // HashTable<label, label, Hash<label> > map(2 * nSize);
// StaticHashTable<label, label, Hash<label> > map(2 * nSize); // StaticHashTable<label, label, Hash<label> > map(2 * nSize);
HashTbl<label, label, Hash<label> > map(2 * nSize); HashTable<label, label, Hash<label> > map(2 * nSize);
Info<< "Constructed map of size: " << nSize Info<< "Constructed map of size: " << nSize
<< " (size " << map.size() << " capacity " << map.capacity() << ") " << " (size " << map.size() << " capacity " << map.capacity() << ") "

14
applications/test/Matrix/MatrixTest.C
View File

@ -24,7 +24,7 @@ License
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "Matrix.H" #include "SquareMatrix.H"
#include "vector.H" #include "vector.H"
using namespace Foam; using namespace Foam;
@ -34,7 +34,7 @@ using namespace Foam;
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
Matrix<scalar> hmm(3, 3); SquareMatrix<scalar> hmm(3);
hmm[0][0] = -3.0; hmm[0][0] = -3.0;
hmm[0][1] = 10.0; hmm[0][1] = 10.0;
@ -46,27 +46,27 @@ int main(int argc, char *argv[])
hmm[2][1] = 6.0; hmm[2][1] = 6.0;
hmm[2][2] = 1.0; hmm[2][2] = 1.0;
Info<< hmm << endl << hmm - 2.0*(-hmm) << endl; //Info<< hmm << endl << hmm - 2.0*(-hmm) << endl;
Info<< max(hmm) << endl; Info<< max(hmm) << endl;
Info<< min(hmm) << endl; Info<< min(hmm) << endl;
Matrix<scalar> hmm2(3, 3, 1.0); SquareMatrix<scalar> hmm2(3, 1.0);
hmm = hmm2; hmm = hmm2;
Info<< hmm << endl; Info<< hmm << endl;
Matrix<scalar> hmm3(Sin); SquareMatrix<scalar> hmm3(Sin);
Info<< hmm3 << endl; Info<< hmm3 << endl;
Matrix<scalar> hmm4; SquareMatrix<scalar> hmm4;
hmm4 = hmm2; hmm4 = hmm2;
Info<< hmm4 << endl; Info<< hmm4 << endl;
Matrix<scalar> hmm5; SquareMatrix<scalar> hmm5;
hmm4 = hmm5; hmm4 = hmm5;
Info<< hmm5 << endl; Info<< hmm5 << endl;

2
applications/test/ODETest/ODETest.C
View File

@ -69,7 +69,7 @@ public:
const scalar x, const scalar x,
const scalarField& y, const scalarField& y,
scalarField& dfdx, scalarField& dfdx,
Matrix<scalar>& dfdy scalarSquareMatrix& dfdy
) const ) const
{ {
dfdx[0] = 0.0; dfdx[0] = 0.0;

2
applications/test/PtrList/PtrListTest.C
View File

@ -57,6 +57,8 @@ public:
Info <<"delete Scalar: " << data_ << endl; Info <<"delete Scalar: " << data_ << endl;
} }
autoPtr<Scalar> clone() const;
friend Ostream& operator<<(Ostream& os, const Scalar& val) friend Ostream& operator<<(Ostream& os, const Scalar& val)
{ {
os << val.data_; os << val.data_;

8
applications/test/UIndirectListTest/UIndirectListTest.C
View File

@ -84,11 +84,11 @@ int main(int argc, char *argv[])
// test List operations // test List operations
List<double> flatList = UIndirectList<double>(completeList, addresses); List<double> flatList(UIndirectList<double>(completeList, addresses));
Info<< "List assigned from UIndirectList: " << flatList << endl; Info<< "List constructed from UIndirectList: " << flatList << endl;
List<double> flatList2(UIndirectList<double>(completeList, addresses)); flatList = UIndirectList<double>(completeList, addresses);
Info<< "List constructed from UIndirectList: " << flatList2 << endl; Info<< "List assigned from UIndirectList: " << flatList << endl;
flatList.append(UIndirectList<double>(completeList, addresses)); flatList.append(UIndirectList<double>(completeList, addresses));
Info<< "List::append(UIndirectList): " << flatList << endl; Info<< "List::append(UIndirectList): " << flatList << endl;

4
applications/test/parallel/parallelTest.C
View File

@ -84,7 +84,7 @@ int main(int argc, char *argv[])
// Sync how many to send // Sync how many to send
labelListList allNTrans(Pstream::nProcs()); labelListList allNTrans(Pstream::nProcs());
allNTrans[Pstream::myProcNo()] = nSend; allNTrans[Pstream::myProcNo()] = nSend;
combineReduce(allNTrans, mapDistribute::listEq()); combineReduce(allNTrans, UPstream::listEq());
// Collect items to be sent // Collect items to be sent
labelListList sendMap(Pstream::nProcs()); labelListList sendMap(Pstream::nProcs());
@ -161,7 +161,7 @@ int main(int argc, char *argv[])
toMaster << data; toMaster << data;
} }
Perr<< "slave receiving from master " Perr<< "slave receiving from master "
<< Pstream::masterNo() << endl; << Pstream::masterNo() << endl;
IPstream fromMaster(Pstream::blocking, Pstream::masterNo()); IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
fromMaster >> data; fromMaster >> data;

14
applications/test/slicedField/slicedFieldTest.C
View File

@ -23,7 +23,7 @@ License
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application Application
volField slicedFieldTest
Description Description
@ -32,16 +32,16 @@ Description
#include "fvCFD.H" #include "fvCFD.H"
#include "SlicedGeometricField.H" #include "SlicedGeometricField.H"
#include "slicedFvPatchFields.H" #include "slicedFvPatchFields.H"
#include "slicedSurfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
#include "setRootCase.H"
# include "setRootCase.H" #include "createTime.H"
#include "createMesh.H"
# include "createTime.H"
# include "createMesh.H"
Info<< "Reading field p\n" << endl; Info<< "Reading field p\n" << endl;
volScalarField p volScalarField p
@ -91,7 +91,7 @@ int main(int argc, char *argv[])
Info<< C << endl; Info<< C << endl;
Info<< (C & U) << endl; Info<< (C & U) << endl;
SlicedGeometricField<vector, fvPatchField, slicedFvPatchField, surfaceMesh> SlicedGeometricField<vector, fvsPatchField, slicedFvsPatchField, surfaceMesh>
Sf Sf
( (
IOobject IOobject
@ -105,7 +105,7 @@ int main(int argc, char *argv[])
mesh.faceAreas() mesh.faceAreas()
); );
Info<< Sf << endl; //Info<< Sf << endl;
return 0; return 0;
} }

6
applications/test/testPointEdgeWave/testPointEdgeWave.C
View File

@ -48,8 +48,6 @@ int main(int argc, char *argv[])
# include "createTime.H" # include "createTime.H"
# include "createPolyMesh.H" # include "createPolyMesh.H"
pointMesh pMesh(mesh);
const polyBoundaryMesh& patches = mesh.boundaryMesh(); const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Get name of patch // Get name of patch
@ -99,7 +97,7 @@ int main(int argc, char *argv[])
PointEdgeWave<pointEdgePoint> wallCalc PointEdgeWave<pointEdgePoint> wallCalc
( (
pMesh, mesh,
wallPoints, wallPoints,
wallInfo, wallInfo,
@ -119,7 +117,7 @@ int main(int argc, char *argv[])
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
), ),
pMesh, pointMesh::New(mesh),
dimensionedScalar("wallDist", dimLength, 0.0) dimensionedScalar("wallDist", dimLength, 0.0)
); );

18
applications/test/volField/volFieldTest.C
View File

@ -34,11 +34,10 @@ Application
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
#include "setRootCase.H"
# include "setRootCase.H" #include "createTime.H"
#include "createMesh.H"
# include "createTime.H"
# include "createMesh.H"
Info<< "Reading field p\n" << endl; Info<< "Reading field p\n" << endl;
volScalarField p volScalarField p
@ -70,12 +69,9 @@ int main(int argc, char *argv[])
mesh mesh
); );
# include "createPhi.H" #include "createPhi.H"
//Info<< transform(dimensionedTensor("I", dimless, 0.1*I), U) << endl; GeometricField<symmTensor, fvPatchField, volMesh> st
GeometricField<sphericalTensor, fvPatchField, volMesh> st
( (
IOobject IOobject
( (
@ -86,8 +82,8 @@ int main(int argc, char *argv[])
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
mesh, mesh,
dimensioned<sphericalTensor>("st", dimless, sphericalTensor::I), dimensioned<symmTensor>("st", dimless, symmTensor::one),
zeroGradientFvPatchSphericalTensorField::typeName zeroGradientFvPatchSymmTensorField::typeName
); );
//Info<< fvc::div(st) << endl; //Info<< fvc::div(st) << endl;

11
applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C
View File

@ -113,7 +113,16 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
return; return;
} }
wordList extensiveVVFNames(IStringStream ("(momentumMean)")()); wordList extensiveVVFNames
(
IStringStream
(
"( \
momentumMean \
fDMean \
)"
)()
);
PtrList<volVectorField> extensiveVVFs(extensiveVVFNames.size()); PtrList<volVectorField> extensiveVVFs(extensiveVVFNames.size());

14
applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C
View File

@ -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; Info<< "Initialising dsmc for Time = " << runTime.timeName() << nl << endl;
dsmcCloud dsmc("dsmc", mesh); dsmcCloud dsmc("dsmc", mesh, dsmcInitialiseDict);
label totalMolecules = dsmc.size(); label totalMolecules = dsmc.size();

1
etc/controlDict
View File

@ -377,6 +377,7 @@ DebugSwitches
displacementLaplacian 0; displacementLaplacian 0;
displacementSBRStress 0; displacementSBRStress 0;
distanceSurface 0; distanceSurface 0;
distribution 0;
downwind 0; downwind 0;
dragModel 0; dragModel 0;
duplicatePoints 0; duplicatePoints 0;

4
src/OpenFOAM/containers/HashTables/HashSet/HashSet.C
View File

@ -105,7 +105,7 @@ template<class Key, class Hash>
void Foam::HashSet<Key, Hash>::operator&=(const HashSet<Key, Hash>& rhs) void Foam::HashSet<Key, Hash>::operator&=(const HashSet<Key, Hash>& rhs)
{ {
// Remove elements not also found in rhs // Remove elements not also found in rhs
for (iterator iter = this->cbegin(); iter != this->cend(); ++iter) for (iterator iter = this->begin(); iter != this->end(); ++iter)
{ {
if (!rhs.found(iter.key())) if (!rhs.found(iter.key()))
{ {
@ -145,8 +145,6 @@ void Foam::HashSet<Key, Hash>::operator-=(const HashSet<Key, Hash>& rhs)
} }
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
/* * * * * * * * * * * * * * * * Global operators * * * * * * * * * * * * * */ /* * * * * * * * * * * * * * * * Global operators * * * * * * * * * * * * * */
template<class Key, class Hash> template<class Key, class Hash>

64
src/OpenFOAM/matrices/Matrix/Matrix.C
View File

@ -235,12 +235,12 @@ void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
template<class Form, class Type> template<class Form, class Type>
const Type& Foam::max(const Matrix<Form, Type>& a) const Type& Foam::max(const Matrix<Form, Type>& a)
{ {
label nm = a.n_*a.m_; label nm = a.n()*a.m();
if (nm) if (nm)
{ {
label curMaxI = 0; label curMaxI = 0;
const Type* v = a.v_[0]; const Type* v = a[0];
for (register label i=1; i<nm; i++) for (register label i=1; i<nm; i++)
{ {
@ -267,12 +267,12 @@ const Type& Foam::max(const Matrix<Form, Type>& a)
template<class Form, class Type> template<class Form, class Type>
const Type& Foam::min(const Matrix<Form, Type>& a) const Type& Foam::min(const Matrix<Form, Type>& a)
{ {
label nm = a.n_*a.m_; label nm = a.n()*a.m();
if (nm) if (nm)
{ {
label curMinI = 0; label curMinI = 0;
const Type* v = a.v_[0]; const Type* v = a[0];
for (register label i=1; i<nm; i++) for (register label i=1; i<nm; i++)
{ {
@ -301,14 +301,14 @@ const Type& Foam::min(const Matrix<Form, Type>& a)
template<class Form, class Type> template<class Form, class Type>
Form Foam::operator-(const Matrix<Form, Type>& a) Form Foam::operator-(const Matrix<Form, Type>& a)
{ {
Form na(a.n_, a.m_); Form na(a.n(), a.m());
if (a.n_ && a.m_) if (a.n() && a.m())
{ {
Type* nav = na.v_[0]; Type* nav = na[0];
const Type* av = a.v_[0]; const Type* av = a[0];
label nm = a.n_*a.m_; label nm = a.n()*a.m();
for (register label i=0; i<nm; i++) for (register label i=0; i<nm; i++)
{ {
nav[i] = -av[i]; nav[i] = -av[i];
@ -322,33 +322,33 @@ Form Foam::operator-(const Matrix<Form, Type>& a)
template<class Form, class Type> template<class Form, class Type>
Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b) Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
{ {
if (a.n_ != b.n_) if (a.n() != b.n())
{ {
FatalErrorIn FatalErrorIn
( (
"Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)" "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
) << "attempted add matrices with different number of rows: " ) << "attempted add matrices with different number of rows: "
<< a.n_ << ", " << b.n_ << a.n() << ", " << b.n()
<< abort(FatalError); << abort(FatalError);
} }
if (a.m_ != b.m_) if (a.m() != b.m())
{ {
FatalErrorIn FatalErrorIn
( (
"Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)" "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
) << "attempted add matrices with different number of columns: " ) << "attempted add matrices with different number of columns: "
<< a.m_ << ", " << b.m_ << a.m() << ", " << b.m()
<< abort(FatalError); << abort(FatalError);
} }
Form ab(a.n_, a.m_); Form ab(a.n(), a.m());
Type* abv = ab.v_[0]; Type* abv = ab[0];
const Type* av = a.v_[0]; const Type* av = a[0];
const Type* bv = b.v_[0]; const Type* bv = b[0];
label nm = a.n_*a.m_; label nm = a.n()*a.m();
for (register label i=0; i<nm; i++) for (register label i=0; i<nm; i++)
{ {
abv[i] = av[i] + bv[i]; abv[i] = av[i] + bv[i];
@ -361,33 +361,33 @@ Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
template<class Form, class Type> template<class Form, class Type>
Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b) Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
{ {
if (a.n_ != b.n_) if (a.n() != b.n())
{ {
FatalErrorIn FatalErrorIn
( (
"Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)" "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
) << "attempted add matrices with different number of rows: " ) << "attempted add matrices with different number of rows: "
<< a.n_ << ", " << b.n_ << a.n() << ", " << b.n()
<< abort(FatalError); << abort(FatalError);
} }
if (a.m_ != b.m_) if (a.m() != b.m())
{ {
FatalErrorIn FatalErrorIn
( (
"Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)" "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
) << "attempted add matrices with different number of columns: " ) << "attempted add matrices with different number of columns: "
<< a.m_ << ", " << b.m_ << a.m() << ", " << b.m()
<< abort(FatalError); << abort(FatalError);
} }
Form ab(a.n_, a.m_); Form ab(a.n(), a.m());
Type* abv = ab.v_[0]; Type* abv = ab[0];
const Type* av = a.v_[0]; const Type* av = a[0];
const Type* bv = b.v_[0]; const Type* bv = b[0];
label nm = a.n_*a.m_; label nm = a.n()*a.m();
for (register label i=0; i<nm; i++) for (register label i=0; i<nm; i++)
{ {
abv[i] = av[i] - bv[i]; abv[i] = av[i] - bv[i];
@ -400,14 +400,14 @@ Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
template<class Form, class Type> template<class Form, class Type>
Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a) Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
{ {
Form sa(a.n_, a.m_); Form sa(a.n(), a.m());
if (a.n_ && a.m_) if (a.n() && a.m())
{ {
Type* sav = sa.v_[0]; Type* sav = sa[0];
const Type* av = a.v_[0]; const Type* av = a[0];
label nm = a.n_*a.m_; label nm = a.n()*a.m();
for (register label i=0; i<nm; i++) for (register label i=0; i<nm; i++)
{ {
sav[i] = s*av[i]; sav[i] = s*av[i];

4
src/OpenFOAM/matrices/Matrix/Matrix.H
View File

@ -161,10 +161,10 @@ public:
// Member operators // Member operators
//- Return subscript-checked element of Matrix. //- Return subscript-checked row of Matrix.
inline Type* operator[](const label); inline Type* operator[](const label);
//- Return subscript-checked element of constant Matrix. //- Return subscript-checked row of constant Matrix.
inline const Type* operator[](const label) const; inline const Type* operator[](const label) const;
//- Assignment operator. Takes linear time. //- Assignment operator. Takes linear time.

3
src/lagrangian/dsmc/Make/files
View File

@ -4,9 +4,6 @@ parcels/derived/dsmcParcel/dsmcParcel.C
/* Cloud base classes */ /* Cloud base classes */
clouds/baseClasses/DsmcBaseCloud/DsmcBaseCloud.C clouds/baseClasses/DsmcBaseCloud/DsmcBaseCloud.C
/* Clouds */
clouds/derived/dsmcCloud/dsmcCloud.C
/* submodels */ /* submodels */
parcels/derived/dsmcParcel/defineDsmcParcel.C parcels/derived/dsmcParcel/defineDsmcParcel.C
parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C

6
src/lagrangian/dsmc/Make/options
View File

@ -1,7 +1,9 @@
EXE_INC = \ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
LIB_LIBS = \ LIB_LIBS = \
-llagrangian \ -llagrangian \
-lfiniteVolume -lfiniteVolume \
-lmeshTools

364
src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
View File

@ -297,7 +297,7 @@ void Foam::DsmcCloud<ParcelType>::collisions()
// Temporary storage for subCells // Temporary storage for subCells
List<DynamicList<label> > subCells(8); List<DynamicList<label> > subCells(8);
scalar deltaT = cachedDeltaT(); scalar deltaT = mesh().time().deltaTValue();
label collisionCandidates = 0; label collisionCandidates = 0;
@ -473,21 +473,92 @@ void Foam::DsmcCloud<ParcelType>::collisions()
template<class ParcelType> template<class ParcelType>
void Foam::DsmcCloud<ParcelType>::resetSurfaceDataFields() void Foam::DsmcCloud<ParcelType>::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<class ParcelType>
void Foam::DsmcCloud<ParcelType>::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<ParcelType>, *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) rhoM *= nParticle_/mesh().cellVolumes();
{ rhoM_.correctBoundaryConditions();
fDBF[p] = vector::zero;
} 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<class ParcelType>
Foam::DsmcCloud<ParcelType>::DsmcCloud Foam::DsmcCloud<ParcelType>::DsmcCloud
( (
const word& cloudName, const word& cloudName,
const volScalarField& T, const fvMesh& mesh,
const volVectorField& U,
bool readFields bool readFields
) )
: :
Cloud<ParcelType>(T.mesh(), cloudName, false), Cloud<ParcelType>(mesh, cloudName, false),
DsmcBaseCloud(), DsmcBaseCloud(),
cloudName_(cloudName), cloudName_(cloudName),
mesh_(T.mesh()), mesh_(mesh),
particleProperties_ particleProperties_
( (
IOobject IOobject
@ -563,37 +633,142 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
( (
IOobject IOobject
( (
this->name() + "q_", "q",
mesh_.time().timeName(), mesh_.time().timeName(),
mesh_, mesh_,
IOobject::NO_READ, IOobject::MUST_READ,
IOobject::NO_WRITE IOobject::AUTO_WRITE
), ),
mesh_, mesh_
dimensionedScalar("zero", dimensionSet(1, 0, -3, 0, 0), 0.0)
), ),
fD_ fD_
( (
IOobject IOobject
( (
this->name() + "fD_", "fD",
mesh_.time().timeName(), mesh_.time().timeName(),
mesh_, mesh_,
IOobject::NO_READ, IOobject::MUST_READ,
IOobject::NO_WRITE IOobject::AUTO_WRITE
), ),
mesh_, mesh_
dimensionedVector ),
rhoN_
(
IOobject
( (
"zero", "rhoN",
dimensionSet(1, -1, -2, 0, 0), mesh_.time().timeName(),
vector::zero 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_(), constProps_(),
rndGen_(label(149382906) + 7183*Pstream::myProcNo()), rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
T_(T), boundaryT_
U_(U), (
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_
( (
BinaryCollisionModel<DsmcCloud<ParcelType> >::New BinaryCollisionModel<DsmcCloud<ParcelType> >::New
@ -641,7 +816,8 @@ template<class ParcelType>
Foam::DsmcCloud<ParcelType>::DsmcCloud Foam::DsmcCloud<ParcelType>::DsmcCloud
( (
const word& cloudName, const word& cloudName,
const fvMesh& mesh const fvMesh& mesh,
const IOdictionary& dsmcInitialiseDict
) )
: :
Cloud<ParcelType>(mesh, cloudName, false), Cloud<ParcelType>(mesh, cloudName, false),
@ -707,15 +883,111 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
vector::zero 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_(), constProps_(),
rndGen_(label(971501) + 1526*Pstream::myProcNo()), rndGen_(label(971501) + 1526*Pstream::myProcNo()),
T_ boundaryT_
( (
volScalarField volScalarField
( (
IOobject IOobject
( (
"T", "boundaryT",
mesh_.time().timeName(), mesh_.time().timeName(),
mesh_, mesh_,
IOobject::NO_READ, IOobject::NO_READ,
@ -725,13 +997,13 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
dimensionedScalar("zero", dimensionSet(0, 0, 0, 1, 0), 0.0) dimensionedScalar("zero", dimensionSet(0, 0, 0, 1, 0), 0.0)
) )
), ),
U_ boundaryU_
( (
volVectorField volVectorField
( (
IOobject IOobject
( (
"U", "boundaryU",
mesh_.time().timeName(), mesh_.time().timeName(),
mesh_, mesh_,
IOobject::NO_READ, IOobject::NO_READ,
@ -754,18 +1026,6 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
buildConstProps(); buildConstProps();
IOdictionary dsmcInitialiseDict
(
IOobject
(
"dsmcInitialiseDict",
mesh_.time().system(),
mesh_,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
initialise(dsmcInitialiseDict); initialise(dsmcInitialiseDict);
} }
@ -782,13 +1042,10 @@ Foam::DsmcCloud<ParcelType>::~DsmcCloud()
template<class ParcelType> template<class ParcelType>
void Foam::DsmcCloud<ParcelType>::evolve() void Foam::DsmcCloud<ParcelType>::evolve()
{ {
// cache the value of deltaT for this timestep
storeDeltaT();
typename ParcelType::trackData td(*this); typename ParcelType::trackData td(*this);
// Reset the surface data collection fields // Reset the data collection fields
resetSurfaceDataFields(); resetFields();
if (debug) if (debug)
{ {
@ -803,6 +1060,9 @@ void Foam::DsmcCloud<ParcelType>::evolve()
// Calculate new velocities via stochastic collisions // Calculate new velocities via stochastic collisions
collisions(); collisions();
// Calculate the volume field data
calculateFields();
} }

143
src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
View File

@ -110,24 +110,41 @@ class DsmcCloud
//- Force density at surface field //- Force density at surface field
volVectorField fD_; 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 //- Parcel constant properties - one for each type
List<typename ParcelType::constantProperties> constProps_; List<typename ParcelType::constantProperties> constProps_;
//- Random number generator //- Random number generator
Random rndGen_; 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 //- boundary velocity
const volScalarField& T_; volVectorField boundaryU_;
//- Velocity
const volVectorField& U_;
// References to the cloud sub-models // References to the cloud sub-models
@ -159,8 +176,11 @@ class DsmcCloud
//- Calculate collisions between molecules //- Calculate collisions between molecules
void collisions(); void collisions();
//- Reset the surface data accumulation field values //- Reset the data accumulation field values to zero
void resetSurfaceDataFields(); void resetFields();
//- Calculate the volume field data
void calculateFields();
//- Disallow default bitwise copy construct //- Disallow default bitwise copy construct
DsmcCloud(const DsmcCloud&); DsmcCloud(const DsmcCloud&);
@ -173,20 +193,21 @@ public:
// Constructors // Constructors
//- Construct given name and mesh, will read Parcels from file //- Construct given name and mesh, will read Parcels and fields from
// file
DsmcCloud DsmcCloud
( (
const word& cloudName, const word& cloudName,
const volScalarField& T, const fvMesh& mesh,
const volVectorField& U,
bool readFields = true bool readFields = true
); );
//- Construct given name and mesh. Used to initialise. //- Construct given name, mesh and initialisation dictionary.
DsmcCloud DsmcCloud
( (
const word& cloudName, const word& cloudName,
const fvMesh& mesh const fvMesh& mesh,
const IOdictionary& dsmcInitialiseDict
); );
@ -242,35 +263,72 @@ public:
//- Return refernce to the random object //- Return refernce to the random object
inline Random& rndGen(); inline Random& rndGen();
//- Store (cache) the current value of deltaT
inline void storeDeltaT();
//- Return the cached value of deltaT // References to the boundary fields for surface data collection
inline scalar cachedDeltaT() const;
//- 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 //- Return non-const number density boundary field reference
inline const volScalarField& q() const; inline volScalarField::GeometricBoundaryField& rhoNBF();
//- Return non-const heat flux at surface field //- Return non-const mass density boundary field reference
inline volScalarField& q(); inline volScalarField::GeometricBoundaryField& rhoMBF();
//- Return force density at surface field //- Return non-const linear kinetic energy density boundary
inline const volVectorField& fD() const; // field reference
inline volScalarField::GeometricBoundaryField& linearKEBF();
//- Return non-const force density at surface field //- Return non-const internal energy density boundary field
inline volVectorField& fD(); // 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 // References to the macroscopic fields
//- Return macroscopic temperature //- Return macroscopic temperature
inline const volScalarField& T() const; inline const volScalarField& boundaryT() const;
//- Return macroscopic velocity //- 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 // Kinetic theory helper functions
@ -385,29 +443,6 @@ public:
void dumpParticlePositions() const; void dumpParticlePositions() const;
// Fields
//- Return the real particle number density field
inline const tmp<volScalarField> rhoN() const;
//- Return the particle mass density field
inline const tmp<volScalarField> rhoM() const;
//- Return the field of number of DSMC particles
inline const tmp<volScalarField> dsmcRhoN() const;
//- Return the momentum density field
inline const tmp<volVectorField> momentum() const;
//- Return the total linear kinetic energy (translational and
// thermal density field
inline const tmp<volScalarField> linearKE() const;
//- Return the internal energy density field
inline const tmp<volScalarField> internalE() const;
//- Return the average internal degrees of freedom field
inline const tmp<volScalarField> iDof() const;
// Cloud evolution functions // Cloud evolution functions

313
src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
View File

@ -126,58 +126,82 @@ inline Foam::Random& Foam::DsmcCloud<ParcelType>::rndGen()
template<class ParcelType> template<class ParcelType>
inline void Foam::DsmcCloud<ParcelType>::storeDeltaT() inline Foam::volScalarField::GeometricBoundaryField&
Foam::DsmcCloud<ParcelType>::qBF()
{ {
cachedDeltaT_ = mesh().time().deltaTValue(); return q_.boundaryField();
} }
template<class ParcelType> template<class ParcelType>
inline Foam::scalar Foam::DsmcCloud<ParcelType>::cachedDeltaT() const inline Foam::volVectorField::GeometricBoundaryField&
Foam::DsmcCloud<ParcelType>::fDBF()
{ {
return cachedDeltaT_; return fD_.boundaryField();
} }
template<class ParcelType> template<class ParcelType>
inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q() const inline Foam::volScalarField::GeometricBoundaryField&
Foam::DsmcCloud<ParcelType>::rhoNBF()
{ {
return q_; return rhoN_.boundaryField();
} }
template<class ParcelType> template<class ParcelType>
inline Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q() inline Foam::volScalarField::GeometricBoundaryField&
Foam::DsmcCloud<ParcelType>::rhoMBF()
{ {
return q_; return rhoM_.boundaryField();
} }
template<class ParcelType> template<class ParcelType>
inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD() const inline Foam::volScalarField::GeometricBoundaryField&
Foam::DsmcCloud<ParcelType>::linearKEBF()
{ {
return fD_; return linearKE_.boundaryField();
} }
template<class ParcelType> template<class ParcelType>
inline Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD() inline Foam::volScalarField::GeometricBoundaryField&
Foam::DsmcCloud<ParcelType>::internalEBF()
{ {
return fD_; return internalE_.boundaryField();
} }
template<class ParcelType> template<class ParcelType>
inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::T() const inline Foam::volScalarField::GeometricBoundaryField&
Foam::DsmcCloud<ParcelType>::iDofBF()
{ {
return T_; return iDof_.boundaryField();
} }
template<class ParcelType> template<class ParcelType>
inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::U() const inline Foam::volVectorField::GeometricBoundaryField&
Foam::DsmcCloud<ParcelType>::momentumBF()
{ {
return U_; return momentum_.boundaryField();
}
template<class ParcelType>
inline const Foam::volScalarField&
Foam::DsmcCloud<ParcelType>::boundaryT() const
{
return boundaryT_;
}
template<class ParcelType>
inline const Foam::volVectorField&
Foam::DsmcCloud<ParcelType>::boundaryU() const
{
return boundaryU_;
} }
@ -381,265 +405,70 @@ Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
template<class ParcelType> template<class ParcelType>
inline const Foam::tmp<Foam::volScalarField> inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q() const
{
return q_;
}
template<class ParcelType>
inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD() const
{
return fD_;
}
template<class ParcelType>
inline const Foam::volScalarField&
Foam::DsmcCloud<ParcelType>::rhoN() const Foam::DsmcCloud<ParcelType>::rhoN() const
{ {
tmp<volScalarField> trhoN return rhoN_;
(
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<ParcelType>, *this, iter)
{
const ParcelType& p = iter();
const label cellI = p.cell();
rhoN[cellI]++;
}
rhoN *= nParticle_/mesh().cellVolumes();
return trhoN;
} }
template<class ParcelType> template<class ParcelType>
inline const Foam::tmp<Foam::volScalarField> inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::rhoM() const
Foam::DsmcCloud<ParcelType>::rhoM() const
{ {
tmp<volScalarField> trhoM return rhoM_;
(
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<ParcelType>, *this, iter)
{
const ParcelType& p = iter();
const label cellI = p.cell();
rhoM[cellI] += constProps(p.typeId()).mass();
}
rhoM *= nParticle_/mesh().cellVolumes();
return trhoM;
} }
template<class ParcelType> template<class ParcelType>
inline const Foam::tmp<Foam::volScalarField> inline const Foam::volScalarField&
Foam::DsmcCloud<ParcelType>::dsmcRhoN() const Foam::DsmcCloud<ParcelType>::dsmcRhoN() const
{ {
tmp<volScalarField> tdsmcRhoN return dsmcRhoN_;
(
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<ParcelType>, *this, iter)
{
const ParcelType& p = iter();
const label cellI = p.cell();
dsmcRhoN[cellI]++;
}
return tdsmcRhoN;
} }
template<class ParcelType> template<class ParcelType>
inline const Foam::tmp<Foam::volVectorField> inline const Foam::volScalarField&
Foam::DsmcCloud<ParcelType>::momentum() const
{
tmp<volVectorField> 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<ParcelType>, *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<class ParcelType>
inline const Foam::tmp<Foam::volScalarField>
Foam::DsmcCloud<ParcelType>::linearKE() const Foam::DsmcCloud<ParcelType>::linearKE() const
{ {
tmp<volScalarField> tlinearKE return linearKE_;
(
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<ParcelType>, *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;
} }
template<class ParcelType> template<class ParcelType>
inline const Foam::tmp<Foam::volScalarField> inline const Foam::volScalarField&
Foam::DsmcCloud<ParcelType>::internalE() const Foam::DsmcCloud<ParcelType>::internalE() const
{ {
tmp<volScalarField> tinternalE return internalE_;
(
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<ParcelType>, *this, iter)
{
const ParcelType& p = iter();
const label cellI = p.cell();
internalE[cellI] += p.Ei();
}
internalE *= nParticle_/mesh().cellVolumes();
return tinternalE;
} }
template<class ParcelType> template<class ParcelType>
inline const Foam::tmp<Foam::volScalarField> inline const Foam::volScalarField&
Foam::DsmcCloud<ParcelType>::iDof() const Foam::DsmcCloud<ParcelType>::iDof() const
{ {
tmp<volScalarField> tiDof return iDof_;
( }
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<ParcelType>, *this, iter) template<class ParcelType>
{ inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::momentum() const
const ParcelType& p = iter(); {
const label cellI = p.cell(); return momentum_;
iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
}
iDof *= nParticle_/mesh().cellVolumes();
return tiDof;
} }

66
src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C
View File

@ -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<dsmcParcel>(cloudName, T, U)
{}
Foam::dsmcCloud::dsmcCloud
(
const word& cloudName,
const fvMesh& mesh
)
:
DsmcCloud<dsmcParcel>(cloudName, mesh)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::dsmcCloud::~dsmcCloud()
{}
// ************************************************************************* //

50
src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H
View File

@ -43,54 +43,8 @@ SourceFiles
namespace Foam namespace Foam
{ {
typedef DsmcCloud<dsmcParcel> dsmcCloud;
/*---------------------------------------------------------------------------*\ }
Class dsmcCloud Declaration
\*---------------------------------------------------------------------------*/
class dsmcCloud
:
public DsmcCloud<dsmcParcel>
{
// 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
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

85
src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
View File

@ -25,6 +25,7 @@ License
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "DsmcParcel.H" #include "DsmcParcel.H"
#include "meshTools.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -43,16 +44,31 @@ bool Foam::DsmcParcel<ParcelType>::move
const polyMesh& mesh = td.cloud().pMesh(); const polyMesh& mesh = td.cloud().pMesh();
const polyBoundaryMesh& pbMesh = mesh.boundaryMesh(); const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
const scalar deltaT = td.cloud().cachedDeltaT(); const scalar deltaT = mesh.time().deltaTValue();
scalar tEnd = (1.0 - p.stepFraction())*deltaT; scalar tEnd = (1.0 - p.stepFraction())*deltaT;
const scalar dtMax = tEnd; 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) 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 // Set the Lagrangian time-step
scalar dt = min(dtMax, tEnd); scalar dt = min(dtMax, tEnd);
dt *= p.trackToFace(p.position() + dt*U_, td); dt *= p.trackToFace(p.position() + dt*Utracking, td);
tEnd -= dt; tEnd -= dt;
@ -113,10 +129,41 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
TrackData& td 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_)); const constantProperties& constProps(td.cloud().constProps(typeId_));
scalar m = constProps.mass(); 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 // pre-interaction energy
scalar preIE = 0.5*m*(U_ & U_) + Ei_; scalar preIE = 0.5*m*(U_ & U_) + Ei_;
@ -132,27 +179,40 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
typeId_ 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 // post-interaction energy
scalar postIE = 0.5*m*(U_ & U_) + Ei_; scalar postIE = 0.5*m*(U_ & U_) + Ei_;
// post-interaction momentum // post-interaction momentum
vector postIMom = m*U_; 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); scalar deltaQ = td.cloud().nParticle()*(preIE - postIE)/(deltaT*fA);
vector deltaFD = td.cloud().nParticle()*(preIMom - postIMom)/(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<ParcelType>::transformProperties
#include "DsmcParcelIO.C" #include "DsmcParcelIO.C"
// ************************************************************************* // // ************************************************************************* //

7
src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C
View File

@ -28,6 +28,7 @@ License
#include "DsmcCloud.H" #include "DsmcCloud.H"
#include "MaxwellianThermal.H" #include "MaxwellianThermal.H"
#include "SpecularReflection.H" #include "SpecularReflection.H"
#include "MixedDiffuseSpecular.H"
namespace Foam namespace Foam
{ {
@ -46,6 +47,12 @@ namespace Foam
DsmcCloud, DsmcCloud,
dsmcParcel dsmcParcel
); );
makeWallInteractionModelType
(
MixedDiffuseSpecular,
DsmcCloud,
dsmcParcel
);
}; };

9
src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
View File

@ -129,7 +129,7 @@ void Foam::FreeStream<CloudType>::inflow()
const polyMesh& mesh(cloud.mesh()); const polyMesh& mesh(cloud.mesh());
const scalar deltaT = cloud.cachedDeltaT(); const scalar deltaT = mesh.time().deltaTValue();
Random& rndGen(cloud.rndGen()); Random& rndGen(cloud.rndGen());
@ -139,12 +139,12 @@ void Foam::FreeStream<CloudType>::inflow()
const volScalarField::GeometricBoundaryField& boundaryT const volScalarField::GeometricBoundaryField& boundaryT
( (
cloud.T().boundaryField() cloud.boundaryT().boundaryField()
); );
const volVectorField::GeometricBoundaryField& boundaryU const volVectorField::GeometricBoundaryField& boundaryU
( (
cloud.U().boundaryField() cloud.boundaryU().boundaryField()
); );
forAll(patches_, p) forAll(patches_, p)
@ -168,7 +168,8 @@ void Foam::FreeStream<CloudType>::inflow()
if (min(boundaryT[patchI]) < SMALL) if (min(boundaryT[patchI]) < SMALL)
{ {
FatalErrorIn ("Foam::FreeStream<CloudType>::inflow()") FatalErrorIn ("Foam::FreeStream<CloudType>::inflow()")
<< "Zero boundary temperature detected, check boundaryT condition." << nl << "Zero boundary temperature detected, check boundaryT "
<< "condition." << nl
<< nl << abort(FatalError); << nl << abort(FatalError);
} }

12
src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
View File

@ -71,10 +71,10 @@ void Foam::MaxwellianThermal<CloudType>::correct
nw /= mag(nw); nw /= mag(nw);
// Normal velocity magnitude // Normal velocity magnitude
scalar magUn = U & nw; scalar U_dot_nw = U & nw;
// Wall tangential velocity (flow direction) // Wall tangential velocity (flow direction)
vector Ut = U - magUn*nw; vector Ut = U - U_dot_nw*nw;
CloudType& cloud(this->owner()); CloudType& cloud(this->owner());
@ -93,9 +93,9 @@ void Foam::MaxwellianThermal<CloudType>::correct
U.z()*(0.8 + 0.2*rndGen.scalar01()) 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 // Wall tangential unit vector
@ -104,7 +104,7 @@ void Foam::MaxwellianThermal<CloudType>::correct
// Other tangential unit vector // Other tangential unit vector
vector tw2 = nw^tw1; vector tw2 = nw^tw1;
scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace]; scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace];
scalar mass = cloud.constProps(typeId).mass(); scalar mass = cloud.constProps(typeId).mass();
@ -118,7 +118,7 @@ void Foam::MaxwellianThermal<CloudType>::correct
- sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw - 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); Ei = cloud.equipartitionInternalEnergy(T, iDof);
} }

140
src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C Normal file
View File

@ -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 <class CloudType>
Foam::MixedDiffuseSpecular<CloudType>::MixedDiffuseSpecular
(
const dictionary& dict,
CloudType& cloud
)
:
WallInteractionModel<CloudType>(dict, cloud, typeName),
diffuseFraction_(readScalar(this->coeffDict().lookup("diffuseFraction")))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template <class CloudType>
Foam::MixedDiffuseSpecular<CloudType>::~MixedDiffuseSpecular()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template <class CloudType>
void Foam::MixedDiffuseSpecular<CloudType>::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;
}
}
}
// ************************************************************************* //

106
src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H Normal file
View File

@ -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 CloudType>
class MixedDiffuseSpecular
:
public WallInteractionModel<CloudType>
{
// 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
// ************************************************************************* //

6
src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
View File

@ -63,11 +63,11 @@ void Foam::SpecularReflection<CloudType>::correct
vector nw = wpp.faceAreas()[wpp.whichFace(faceId)]; vector nw = wpp.faceAreas()[wpp.whichFace(faceId)];
nw /= mag(nw); 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;
} }
} }

108
src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C
View File

@ -25,27 +25,49 @@ License
\*----------------------------------------------------------------------------*/ \*----------------------------------------------------------------------------*/
#include "distribution.H" #include "distribution.H"
#include "OFstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam namespace Foam
{ {
defineTypeNameAndDebug(distribution, 0);
}
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
void Foam::distribution::write
(
const fileName& file,
const List<Pair<scalar> >& pairs
)
{
OFstream os(file);
forAll(pairs, i)
{
os << pairs[i].first() << ' ' << pairs[i].second() << nl;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
distribution::distribution() Foam::distribution::distribution()
: :
Map<label>(), Map<label>(),
binWidth_(1) binWidth_(1)
{} {}
distribution::distribution(const scalar binWidth) Foam::distribution::distribution(const scalar binWidth)
: :
Map<label>(), Map<label>(),
binWidth_(binWidth) binWidth_(binWidth)
{} {}
distribution::distribution(const distribution& d) Foam::distribution::distribution(const distribution& d)
: :
Map<label>(static_cast< Map<label> >(d)), Map<label>(static_cast< Map<label> >(d)),
binWidth_(d.binWidth()) binWidth_(d.binWidth())
@ -54,13 +76,13 @@ distribution::distribution(const distribution& d)
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
distribution::~distribution() Foam::distribution::~distribution()
{} {}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
label distribution::totalEntries() const Foam::label Foam::distribution::totalEntries() const
{ {
label sumOfEntries = 0; label sumOfEntries = 0;
@ -88,7 +110,7 @@ label distribution::totalEntries() const
} }
scalar distribution::approxTotalEntries() const Foam::scalar Foam::distribution::approxTotalEntries() const
{ {
scalar sumOfEntries = 0; scalar sumOfEntries = 0;
@ -101,7 +123,7 @@ scalar distribution::approxTotalEntries() const
} }
scalar distribution::mean() const Foam::scalar Foam::distribution::mean() const
{ {
scalar runningSum = 0; scalar runningSum = 0;
@ -124,7 +146,7 @@ scalar distribution::mean() const
} }
scalar distribution::median() Foam::scalar Foam::distribution::median()
{ {
// From: // From:
// http://mathworld.wolfram.com/StatisticalMedian.html // http://mathworld.wolfram.com/StatisticalMedian.html
@ -188,7 +210,7 @@ scalar distribution::median()
} }
void distribution::add(const scalar valueToAdd) void Foam::distribution::add(const scalar valueToAdd)
{ {
iterator iter(this->begin()); iterator iter(this->begin());
@ -218,13 +240,13 @@ void distribution::add(const scalar valueToAdd)
} }
void distribution::add(const label valueToAdd) void Foam::distribution::add(const label valueToAdd)
{ {
add(scalar(valueToAdd)); add(scalar(valueToAdd));
} }
void distribution::insertMissingKeys() void Foam::distribution::insertMissingKeys()
{ {
iterator iter(this->begin()); iterator iter(this->begin());
@ -247,7 +269,7 @@ void distribution::insertMissingKeys()
} }
List< Pair<scalar> > distribution::normalised() Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalised()
{ {
scalar totEnt = approxTotalEntries(); scalar totEnt = approxTotalEntries();
@ -268,17 +290,25 @@ List< Pair<scalar> > distribution::normalised()
normDist[k].second() = scalar((*this)[key])/totEnt/binWidth_; normDist[k].second() = scalar((*this)[key])/totEnt/binWidth_;
} }
if (debug)
{
Info<< "totEnt: " << totEnt << endl;
}
return normDist; return normDist;
} }
List< Pair<scalar> > distribution::normalisedMinusMean() Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalisedMinusMean()
{ {
return normalisedShifted(mean()); return normalisedShifted(mean());
} }
List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue) Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalisedShifted
(
scalar shiftValue
)
{ {
List<Pair<scalar> > oldDist(normalised()); List<Pair<scalar> > oldDist(normalised());
@ -301,20 +331,23 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
label newKey = lowestNewKey; label newKey = lowestNewKey;
// Info << shiftValue if (debug)
// << nl << lowestOldBin {
// << nl << lowestNewKey Info<< shiftValue
// << nl << interpolationStartDirection << nl << lowestOldBin
// << endl; << nl << lowestNewKey
<< nl << interpolationStartDirection
<< endl;
// scalar checkNormalisation = 0; scalar checkNormalisation = 0;
// forAll (oldDist, oD) forAll (oldDist, oD)
// { {
// checkNormalisation += oldDist[oD].second()*binWidth_; checkNormalisation += oldDist[oD].second()*binWidth_;
// } }
// Info << "Initial normalisation = " << checkNormalisation << endl; Info<< "Initial normalisation = " << checkNormalisation << endl;
}
forAll(oldDist,u) forAll(oldDist,u)
{ {
@ -368,20 +401,23 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
newKey++; newKey++;
} }
// checkNormalisation = 0; if (debug)
{
scalar checkNormalisation = 0;
// forAll (newDist, nD) forAll (newDist, nD)
// { {
// checkNormalisation += newDist[nD].second()*binWidth_; checkNormalisation += newDist[nD].second()*binWidth_;
// } }
// Info << "Shifted normalisation = " << checkNormalisation << endl; Info<< "Shifted normalisation = " << checkNormalisation << endl;
}
return newDist; return newDist;
} }
List<Pair<scalar> > distribution::raw() Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::raw()
{ {
insertMissingKeys(); insertMissingKeys();
@ -406,7 +442,7 @@ List<Pair<scalar> > distribution::raw()
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void distribution::operator=(const distribution& rhs) void Foam::distribution::operator=(const distribution& rhs)
{ {
// Check for assignment to self // Check for assignment to self
if (this == &rhs) if (this == &rhs)
@ -424,7 +460,7 @@ void distribution::operator=(const distribution& rhs)
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const distribution& d) Foam::Ostream& Foam::operator<<(Ostream& os, const distribution& d)
{ {
os << d.binWidth_ os << d.binWidth_
<< static_cast<const Map<label>&>(d); << static_cast<const Map<label>&>(d);
@ -440,8 +476,4 @@ Ostream& operator<<(Ostream& os, const distribution& d)
} }
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* // // ************************************************************************* //

21
src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.H
View File

@ -26,6 +26,8 @@ Class
Foam::distribution Foam::distribution
Description Description
Accumulating histogram of values. Specified bin resolution
automatic generation of bins.
SourceFiles SourceFiles
distributionI.H distributionI.H
@ -59,6 +61,21 @@ class distribution
public: public:
//- Runtime type information
TypeName("distribution");
// Static functions
//- write to file
static void write
(
const fileName& file,
const List<Pair<scalar> >& pairs
);
// Constructors // Constructors
//- Construct null //- Construct null
@ -73,7 +90,7 @@ public:
// Destructor // Destructor
~distribution(); virtual ~distribution();
// Member Functions // Member Functions
@ -97,7 +114,7 @@ public:
List<Pair<scalar> > normalisedMinusMean(); List<Pair<scalar> > normalisedMinusMean();
List<Pair<scalar> > normalisedShifted(const scalar shiftValue); List<Pair<scalar> > normalisedShifted(scalar shiftValue);
List<Pair<scalar> > raw(); List<Pair<scalar> > raw();

46
src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
View File

@ -110,6 +110,7 @@ void Foam::dsmcFields::write()
word linearKEMeanName = "linearKEMean"; word linearKEMeanName = "linearKEMean";
word internalEMeanName = "internalEMean"; word internalEMeanName = "internalEMean";
word iDofMeanName = "iDofMean"; word iDofMeanName = "iDofMean";
word fDMeanName = "fDMean";
const volScalarField& rhoNMean = obr_.lookupObject<volScalarField> const volScalarField& rhoNMean = obr_.lookupObject<volScalarField>
( (
@ -141,6 +142,11 @@ void Foam::dsmcFields::write()
iDofMeanName iDofMeanName
); );
volVectorField fDMean = obr_.lookupObject<volVectorField>
(
fDMeanName
);
if (min(mag(rhoNMean)).value() > VSMALL) if (min(mag(rhoNMean)).value() > VSMALL)
{ {
Info<< "Calculating dsmcFields." << endl; Info<< "Calculating dsmcFields." << endl;
@ -168,8 +174,9 @@ void Foam::dsmcFields::write()
obr_, obr_,
IOobject::NO_READ IOobject::NO_READ
), ),
2.0/(3.0*physicoChemical::k.value()*rhoNMean) 2.0/(3.0*physicoChemical::k.value()*rhoNMean)
*(linearKEMean - 0.5*rhoMMean*(UMean & UMean)) *(linearKEMean - 0.5*rhoMMean*(UMean & UMean))
); );
Info<< " Calculating internalT field." << endl; Info<< " Calculating internalT field." << endl;
@ -182,7 +189,7 @@ void Foam::dsmcFields::write()
obr_, obr_,
IOobject::NO_READ IOobject::NO_READ
), ),
2.0/(physicoChemical::k.value()*iDofMean)*internalEMean (2.0/physicoChemical::k.value())*(internalEMean/iDofMean)
); );
Info<< " Calculating overallT field." << endl; Info<< " Calculating overallT field." << endl;
@ -196,9 +203,36 @@ void Foam::dsmcFields::write()
IOobject::NO_READ IOobject::NO_READ
), ),
2.0/(physicoChemical::k.value()*(3.0*rhoNMean + iDofMean)) 2.0/(physicoChemical::k.value()*(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
),
physicoChemical::k.value()*rhoNMean*translationalT
);
const fvMesh& mesh = fDMean.mesh();
forAll(mesh.boundaryMesh(), i)
{
const polyPatch& patch = mesh.boundaryMesh()[i];
if (isA<wallPolyPatch>(patch))
{
p.boundaryField()[i] =
fDMean.boundaryField()[i]
& (patch.faceAreas()/mag(patch.faceAreas()));
}
}
Info<< " mag(UMean) max/min : " Info<< " mag(UMean) max/min : "
<< max(mag(UMean)).value() << " " << max(mag(UMean)).value() << " "
<< min(mag(UMean)).value() << endl; << min(mag(UMean)).value() << endl;
@ -215,6 +249,10 @@ void Foam::dsmcFields::write()
<< max(overallT).value() << " " << max(overallT).value() << " "
<< min(overallT).value() << endl; << min(overallT).value() << endl;
Info<< " p max/min : "
<< max(p).value() << " "
<< min(p).value() << endl;
UMean.write(); UMean.write();
translationalT.write(); translationalT.write();
@ -223,6 +261,8 @@ void Foam::dsmcFields::write()
overallT.write(); overallT.write();
p.write();
Info<< "dsmcFields written." << nl << endl; Info<< "dsmcFields written." << nl << endl;
} }
else else

3
tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/iDof
View File

@ -32,7 +32,8 @@ boundaryField
walls walls
{ {
type zeroGradient; type calculated;
value uniform 0;
} }
} }

3
tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/internalE
View File

@ -32,7 +32,8 @@ boundaryField
walls walls
{ {
type zeroGradient; type calculated;
value uniform 0;
} }
} }

3
tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/linearKE
View File

@ -32,7 +32,8 @@ boundaryField
walls walls
{ {
type zeroGradient; type calculated;
value uniform 0;
} }
} }

3
tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/momentum
View File

@ -32,7 +32,8 @@ boundaryField
walls walls
{ {
type zeroGradient; type calculated;
value uniform (0 0 0);
} }
} }

3
tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/rhoM
View File

@ -32,7 +32,8 @@ boundaryField
walls walls
{ {
type zeroGradient; type calculated;
value uniform 0;
} }
} }

3
tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/rhoN
View File

@ -32,7 +32,8 @@ boundaryField
walls walls
{ {
type zeroGradient; type calculated;
value uniform 0;
} }
} }

3
tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/iDof
View File

@ -27,7 +27,8 @@ boundaryField
obstacle obstacle
{ {
type zeroGradient; type calculated;
value uniform 0;
} }
periodic periodic

3
tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/internalE
View File

@ -27,7 +27,8 @@ boundaryField
obstacle obstacle
{ {
type zeroGradient; type calculated;
value uniform 0;
} }
periodic periodic

3
tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/linearKE
View File

@ -27,7 +27,8 @@ boundaryField
obstacle obstacle
{ {
type zeroGradient; type calculated;
value uniform 0;
} }
periodic periodic

3
tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/momentum
View File

@ -27,7 +27,8 @@ boundaryField
obstacle obstacle
{ {
type zeroGradient; type calculated;
value uniform (0 0 0);
} }
periodic periodic

3
tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/rhoM
View File

@ -27,7 +27,8 @@ boundaryField
obstacle obstacle
{ {
type zeroGradient; type calculated;
value uniform 0;
} }
periodic periodic

3
tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/rhoN
View File

@ -27,7 +27,8 @@ boundaryField
obstacle obstacle
{ {
type zeroGradient; type calculated;
value uniform 0;
} }
periodic periodic