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 branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
mattijs
2013-02-26 11:26:20 +00:00
parent d823a0be75 0a5508e813
commit 1d20c225c3
40 changed files with 1342 additions and 1464 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/finiteVolume/Make/files
View File

@ -340,8 +340,7 @@ $(gradSchemes)/gaussGrad/gaussGrads.C
$(gradSchemes)/leastSquaresGrad/leastSquaresVectors.C
$(gradSchemes)/leastSquaresGrad/leastSquaresGrads.C
$(gradSchemes)/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C
$(gradSchemes)/extendedLeastSquaresGrad/extendedLeastSquaresGrads.C
$(gradSchemes)/LeastSquaresGrad/LeastSquaresGrads.C
$(gradSchemes)/fourthGrad/fourthGrads.C
limitedGradSchemes = $(gradSchemes)/limitedGradSchemes

4
src/finiteVolume/fields/fvPatchFields/derived/turbulentInlet/turbulentInletFvPatchField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -39,7 +39,7 @@ Description
\vartable
x_p | patch values
x_{ref} | refernce patch values
x_{ref} | reference patch values
n | time level
\alpha | fraction of new random component added to previous time value
C_{RMS} | RMS coefficient

115
src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.C → src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresGrad.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,18 +23,16 @@ License
\*---------------------------------------------------------------------------*/
#include "extendedLeastSquaresGrad.H"
#include "extendedLeastSquaresVectors.H"
#include "LeastSquaresGrad.H"
#include "LeastSquaresVectors.H"
#include "gaussGrad.H"
#include "fvMesh.H"
#include "volMesh.H"
#include "surfaceMesh.H"
#include "GeometricField.H"
#include "zeroGradientFvPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
template<class Type, class Stencil>
Foam::tmp
<
Foam::GeometricField
@ -44,15 +42,21 @@ Foam::tmp
Foam::volMesh
>
>
Foam::fv::extendedLeastSquaresGrad<Type>::calcGrad
Foam::fv::LeastSquaresGrad<Type, Stencil>::calcGrad
(
const GeometricField<Type, fvPatchField, volMesh>& vsf,
const GeometricField<Type, fvPatchField, volMesh>& vtf,
const word& name
) const
{
typedef typename outerProduct<vector, Type>::type GradType;
const fvMesh& mesh = vsf.mesh();
const fvMesh& mesh = vtf.mesh();
// Get reference to least square vectors
const LeastSquaresVectors<Stencil>& lsv = LeastSquaresVectors<Stencil>::New
(
mesh
);
tmp<GeometricField<GradType, fvPatchField, volMesh> > tlsGrad
(
@ -61,7 +65,7 @@ Foam::fv::extendedLeastSquaresGrad<Type>::calcGrad
IOobject
(
name,
vsf.instance(),
vtf.instance(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
@ -70,87 +74,64 @@ Foam::fv::extendedLeastSquaresGrad<Type>::calcGrad
dimensioned<GradType>
(
"zero",
vsf.dimensions()/dimLength,
vtf.dimensions()/dimLength,
pTraits<GradType>::zero
),
zeroGradientFvPatchField<GradType>::typeName
)
);
GeometricField<GradType, fvPatchField, volMesh>& lsGrad = tlsGrad();
Field<GradType>& lsGradIf = lsGrad;
// Get reference to least square vectors
const extendedLeastSquaresVectors& lsv = extendedLeastSquaresVectors::New
(
mesh,
minDet_
);
const extendedCentredCellToCellStencil& stencil = lsv.stencil();
const List<List<label> >& stencilAddr = stencil.stencil();
const List<List<vector> >& lsvs = lsv.vectors();
const surfaceVectorField& ownLs = lsv.pVectors();
const surfaceVectorField& neiLs = lsv.nVectors();
// Construct flat version of vtf
// including all values referred to by the stencil
List<Type> flatVtf(stencil.map().constructSize(), pTraits<Type>::zero);
const labelUList& owner = mesh.owner();
const labelUList& neighbour = mesh.neighbour();
forAll(owner, facei)
// Insert internal values
forAll(vtf, celli)
{
label own = owner[facei];
label nei = neighbour[facei];
Type deltaVsf = vsf[nei] - vsf[own];
lsGrad[own] += ownLs[facei]*deltaVsf;
lsGrad[nei] -= neiLs[facei]*deltaVsf;
flatVtf[celli] = vtf[celli];
}
// Boundary faces
forAll(vsf.boundaryField(), patchi)
// Insert boundary values
forAll(vtf.boundaryField(), patchi)
{
const fvsPatchVectorField& patchOwnLs = ownLs.boundaryField()[patchi];
const fvPatchField<Type>& ptf = vtf.boundaryField()[patchi];
const labelUList& faceCells =
lsGrad.boundaryField()[patchi].patch().faceCells();
label nCompact =
ptf.patch().start()
- mesh.nInternalFaces()
+ mesh.nCells();
if (vsf.boundaryField()[patchi].coupled())
forAll(ptf, i)
{
const Field<Type> neiVsf
(
vsf.boundaryField()[patchi].patchNeighbourField()
);
forAll(neiVsf, patchFaceI)
{
lsGrad[faceCells[patchFaceI]] +=
patchOwnLs[patchFaceI]
*(neiVsf[patchFaceI] - vsf[faceCells[patchFaceI]]);
}
}
else
{
const fvPatchField<Type>& patchVsf = vsf.boundaryField()[patchi];
forAll(patchVsf, patchFaceI)
{
lsGrad[faceCells[patchFaceI]] +=
patchOwnLs[patchFaceI]
*(patchVsf[patchFaceI] - vsf[faceCells[patchFaceI]]);
}
flatVtf[nCompact++] = ptf[i];
}
}
// Do all swapping to complete flatVtf
stencil.map().distribute(flatVtf);
const List<labelPair>& additionalCells = lsv.additionalCells();
const vectorField& additionalVectors = lsv.additionalVectors();
forAll(additionalCells, i)
// Accumulate the cell-centred gradient from the
// weighted least-squares vectors and the flattened field values
forAll(stencilAddr, celli)
{
lsGrad[additionalCells[i][0]] +=
additionalVectors[i]
*(vsf[additionalCells[i][1]] - vsf[additionalCells[i][0]]);
const labelList& compactCells = stencilAddr[celli];
const List<vector>& lsvc = lsvs[celli];
forAll(compactCells, i)
{
lsGradIf[celli] += lsvc[i]*flatVtf[compactCells[i]];
}
}
// Correct the boundary conditions
lsGrad.correctBoundaryConditions();
gaussGrad<Type>::correctBoundaryConditions(vsf, lsGrad);
gaussGrad<Type>::correctBoundaryConditions(vtf, lsGrad);
return tlsGrad;
}

174
src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresGrad.H Normal file
View File

@ -0,0 +1,174 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::fv::LeastSquaresGrad
Description
Gradient calculated using weighted least-squares on an arbitrary stencil.
The stencil type is provided via a template argument and any cell-based
stencil is supported:
\table
Stencil | Connections | Scheme name
centredCFCCellToCellStencil | cell-face-cell | Not Instantiated
centredCPCCellToCellStencil | cell-point-cell | pointCellsLeastSquares
centredCECCellToCellStencil | cell-edge-cell | edgeCellsLeastSquares
\endtable
The first of these is not instantiated by default as the standard
leastSquaresGrad is equivalent and more efficient.
\heading Usage
Example of the gradient specification:
\verbatim
gradSchemes
{
default pointCellsLeastSquares;
}
\endverbatim
See Also
Foam::fv::LeastSquaresVectors
Foam::fv::leastSquaresGrad
SourceFiles
LeastSquaresGrad.C
LeastSquaresVectors.H
LeastSquaresVectors.C
\*---------------------------------------------------------------------------*/
#ifndef LeastSquaresGrad_H
#define LeastSquaresGrad_H
#include "gradScheme.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fv
{
/*---------------------------------------------------------------------------*\
Class LeastSquaresGrad Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class Stencil>
class LeastSquaresGrad
:
public fv::gradScheme<Type>
{
// Private Data
//- Minimum determinant criterion to choose extra cells
scalar minDet_;
// Private Member Functions
//- Disallow default bitwise copy construct
LeastSquaresGrad(const LeastSquaresGrad&);
//- Disallow default bitwise assignment
void operator=(const LeastSquaresGrad&);
public:
//- Runtime type information
TypeName("LeastSquares");
// Constructors
//- Construct from Istream
LeastSquaresGrad(const fvMesh& mesh, Istream& schemeData)
:
gradScheme<Type>(mesh)
{}
// Member Functions
//- Return the gradient of the given field to the gradScheme::grad
// for optional caching
virtual tmp
<
GeometricField
<typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
> calcGrad
(
const GeometricField<Type, fvPatchField, volMesh>& vsf,
const word& name
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Add the patch constructor functions to the hash tables
#define makeLeastSquaresGradTypeScheme(SS, STENCIL, TYPE) \
\
typedef LeastSquaresGrad<TYPE, STENCIL> LeastSquaresGrad##TYPE##STENCIL##_; \
defineTemplateTypeNameAndDebugWithName \
(LeastSquaresGrad##TYPE##STENCIL##_, #SS, 0); \
\
gradScheme<TYPE>::addIstreamConstructorToTable \
<LeastSquaresGrad<TYPE, STENCIL> > \
add##SS##STENCIL##TYPE##IstreamConstructorToTable_;
#define makeLeastSquaresGradScheme(SS, STENCIL) \
\
typedef LeastSquaresVectors<STENCIL> LeastSquaresVectors##STENCIL##_; \
defineTemplateTypeNameAndDebugWithName \
(LeastSquaresVectors##STENCIL##_, #SS, 0); \
\
makeLeastSquaresGradTypeScheme(SS,STENCIL,scalar) \
makeLeastSquaresGradTypeScheme(SS,STENCIL,vector)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "LeastSquaresGrad.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

57
src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresGrads.C Normal file
View File

@ -0,0 +1,57 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "LeastSquaresGrad.H"
//#include "centredCFCCellToCellStencilObject.H"
#include "centredCPCCellToCellStencilObject.H"
#include "centredCECCellToCellStencilObject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace fv
{
// makeLeastSquaresGradScheme
// (
// faceCellsLeastSquares,
// centredCFCCellToCellStencilObject
// )
makeLeastSquaresGradScheme
(
pointCellsLeastSquares,
centredCPCCellToCellStencilObject
)
makeLeastSquaresGradScheme
(
edgeCellsLeastSquares,
centredCECCellToCellStencilObject
)
}
}
// ************************************************************************* //

118
src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresVectors.C Normal file
View File

@ -0,0 +1,118 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "LeastSquaresVectors.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Stencil>
Foam::fv::LeastSquaresVectors<Stencil>::LeastSquaresVectors
(
const fvMesh& mesh
)
:
MeshObject<fvMesh, LeastSquaresVectors>(mesh),
vectors_(mesh.nCells())
{
calcLeastSquaresVectors();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Stencil>
Foam::fv::LeastSquaresVectors<Stencil>::~LeastSquaresVectors()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Stencil>
void Foam::fv::LeastSquaresVectors<Stencil>::calcLeastSquaresVectors()
{
if (debug)
{
Info<< "LeastSquaresVectors::calcLeastSquaresVectors() :"
<< "Calculating least square gradient vectors"
<< endl;
}
const fvMesh& mesh = this->mesh_;
const extendedCentredCellToCellStencil& stencil = this->stencil();
stencil.collectData(mesh.C(), vectors_);
// Create the base form of the dd-tensor
// including components for the "empty" directions
symmTensor dd0(sqr((Vector<label>::one - mesh.geometricD())/2));
forAll (vectors_, i)
{
List<vector>& lsvi = vectors_[i];
symmTensor dd(dd0);
// The current cell is 0 in the stencil
// Calculate the deltas and sum the weighted dd
for (label j=1; j<lsvi.size(); j++)
{
lsvi[j] = lsvi[j] - lsvi[0];
scalar magSqrLsvi = magSqr(lsvi[j]);
dd += sqr(lsvi[j])/magSqrLsvi;
lsvi[j] /= magSqrLsvi;
}
// Invert dd
dd = inv(dd);
// Remove the components corresponding to the empty directions
dd -= dd0;
// Finalize the gradient weighting vectors
lsvi[0] = vector::zero;
for (label j=1; j<lsvi.size(); j++)
{
lsvi[j] = dd & lsvi[j];
lsvi[0] -= lsvi[j];
}
}
if (debug)
{
Info<< "LeastSquaresVectors::calcLeastSquaresVectors() :"
<< "Finished calculating least square gradient vectors"
<< endl;
}
}
template<class Stencil>
bool Foam::fv::LeastSquaresVectors<Stencil>::movePoints()
{
calcLeastSquaresVectors();
return true;
}
// ************************************************************************* //

85
src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.H → src/finiteVolume/finiteVolume/gradSchemes/LeastSquaresGrad/LeastSquaresVectors.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,103 +22,110 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::extendedLeastSquaresVectors
Foam::fv::LeastSquaresVectors
Description
Extended molecule least-squares gradient scheme vectors
Least-squares gradient scheme vectors
See Also
Foam::fv::LeastSquaresGrad
SourceFiles
extendedLeastSquaresVectors.C
LeastSquaresVectors.C
\*---------------------------------------------------------------------------*/
#ifndef extendedLeastSquaresVectors_H
#define extendedLeastSquaresVectors_H
#ifndef LeastSquaresVectors_H
#define LeastSquaresVectors_H
#include "extendedCentredCellToCellStencil.H"
#include "MeshObject.H"
#include "fvMesh.H"
#include "surfaceFieldsFwd.H"
#include "labelPair.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fv
{
/*---------------------------------------------------------------------------*\
Class extendedLeastSquaresVectors Declaration
Class LeastSquaresVectors Declaration
\*---------------------------------------------------------------------------*/
class extendedLeastSquaresVectors
template<class Stencil>
class LeastSquaresVectors
:
public MeshObject<fvMesh, extendedLeastSquaresVectors>
public MeshObject<fvMesh, LeastSquaresVectors<Stencil> >
{
// Private data
//- Minimum determinant criterion to choose extra cells
scalar minDet_;
//- Least-squares gradient vectors
mutable surfaceVectorField* pVectorsPtr_;
mutable surfaceVectorField* nVectorsPtr_;
mutable List<labelPair>* additionalCellsPtr_;
mutable vectorField* additionalVectorsPtr_;
List<List<vector> > vectors_;
// Private Member Functions
//- Construct Least-squares gradient vectors
void makeLeastSquaresVectors() const;
//- Calculate Least-squares gradient vectors
void calcLeastSquaresVectors();
public:
// Declare name of the class and its debug switch
TypeName("extendedLeastSquaresVectors");
TypeName("LeastSquaresVectors");
// Constructors
//- Construct given an fvMesh and the minimum determinant criterion
extendedLeastSquaresVectors
LeastSquaresVectors
(
const fvMesh&,
const scalar minDet
const fvMesh&
);
//- Destructor
virtual ~extendedLeastSquaresVectors();
virtual ~LeastSquaresVectors();
// Member functions
//- Return reference to owner least square vectors
const surfaceVectorField& pVectors() const;
//- Return reference to neighbour least square vectors
const surfaceVectorField& nVectors() const;
//- Return reference to additional cells
const List<labelPair>& additionalCells() const;
//- Return reference to additional least square vectors
const vectorField& additionalVectors() const;
//- Return reference to the stencil
const extendedCentredCellToCellStencil& stencil() const
{
return Stencil::New(this->mesh_);
}
//- Return reference to the least square vectors
const List<List<vector> >& vectors() const
{
return vectors_;
}
//- Delete the least square vectors when the mesh moves
virtual bool movePoints();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "LeastSquaresVectors.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

136
src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrad.H
View File

@ -1,136 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::fv::extendedLeastSquaresGrad
Description
Second-order gradient scheme using least-squares.
SourceFiles
extendedLeastSquaresGrad.C
\*---------------------------------------------------------------------------*/
#ifndef extendedLeastSquaresGrad_H
#define extendedLeastSquaresGrad_H
#include "gradScheme.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fv
{
/*---------------------------------------------------------------------------*\
Class extendedLeastSquaresGrad Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class extendedLeastSquaresGrad
:
public fv::gradScheme<Type>
{
// Private Data
//- Minimum determinant criterion to choose extra cells
scalar minDet_;
// Private Member Functions
//- Disallow default bitwise copy construct
extendedLeastSquaresGrad(const extendedLeastSquaresGrad&);
//- Disallow default bitwise assignment
void operator=(const extendedLeastSquaresGrad&);
public:
//- Runtime type information
TypeName("extendedLeastSquares");
// Constructors
//- Construct from Istream
extendedLeastSquaresGrad(const fvMesh& mesh, Istream& schemeData)
:
gradScheme<Type>(mesh),
minDet_(readScalar(schemeData))
{
if (minDet_ < 0) //-for facearea weighted: || minDet_ > 8)
{
FatalIOErrorIn
(
"extendedLeastSquaresGrad"
"(const fvMesh&, Istream& schemeData)",
schemeData
) << "Minimum determinant = " << minDet_
<< " should be >= 0" // and <= 8"
<< exit(FatalIOError);
}
}
// Member Functions
//- Return the gradient of the given field to the gradScheme::grad
// for optional caching
virtual tmp
<
GeometricField
<typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
> calcGrad
(
const GeometricField<Type, fvPatchField, volMesh>& vsf,
const word& name
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "extendedLeastSquaresGrad.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

453
src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C
View File

@ -1,453 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "extendedLeastSquaresVectors.H"
#include "surfaceFields.H"
#include "volFields.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(extendedLeastSquaresVectors, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::extendedLeastSquaresVectors::extendedLeastSquaresVectors
(
const fvMesh& mesh,
const scalar minDet
)
:
MeshObject<fvMesh, extendedLeastSquaresVectors>(mesh),
minDet_(minDet),
pVectorsPtr_(NULL),
nVectorsPtr_(NULL),
additionalCellsPtr_(NULL),
additionalVectorsPtr_(NULL)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::extendedLeastSquaresVectors::~extendedLeastSquaresVectors()
{
deleteDemandDrivenData(pVectorsPtr_);
deleteDemandDrivenData(nVectorsPtr_);
deleteDemandDrivenData(additionalCellsPtr_);
deleteDemandDrivenData(additionalVectorsPtr_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::extendedLeastSquaresVectors::makeLeastSquaresVectors() const
{
if (debug)
{
Info<< "extendedLeastSquaresVectors::makeLeastSquaresVectors() :"
<< "Constructing least square gradient vectors"
<< endl;
}
const fvMesh& mesh = mesh_;
pVectorsPtr_ = new surfaceVectorField
(
IOobject
(
"LeastSquaresP",
mesh_.pointsInstance(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh_,
dimensionedVector("zero", dimless/dimLength, vector::zero)
);
surfaceVectorField& lsP = *pVectorsPtr_;
nVectorsPtr_ = new surfaceVectorField
(
IOobject
(
"LeastSquaresN",
mesh_.pointsInstance(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh_,
dimensionedVector("zero", dimless/dimLength, vector::zero)
);
surfaceVectorField& lsN = *nVectorsPtr_;
// Set local references to mesh data
const labelUList& owner = mesh_.owner();
const labelUList& neighbour = mesh_.neighbour();
// Determine number of dimensions and (for 2D) missing dimension
label nDims = 0;
label twoD = -1;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (mesh.geometricD()[dir] == 1)
{
nDims++;
}
else
{
twoD = dir;
}
}
if (nDims == 1)
{
FatalErrorIn
(
"extendedLeastSquaresVectors::makeLeastSquaresVectors() const"
) << "Found a mesh with only one geometric dimension : "
<< mesh.geometricD()
<< exit(FatalError);
}
else if (nDims == 2)
{
Info<< "extendedLeastSquares : detected " << nDims
<< " valid directions. Missing direction " << twoD << endl;
}
const volVectorField& C = mesh.C();
// Set up temporary storage for the dd tensor (before inversion)
symmTensorField dd(mesh_.nCells(), symmTensor::zero);
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
vector d = C[nei] - C[own];
const symmTensor wdd(1.0/magSqr(d)*sqr(d));
dd[own] += wdd;
dd[nei] += wdd;
}
// Visit the boundaries. Coupled boundaries are taken into account
// in the construction of d vectors.
surfaceVectorField::GeometricBoundaryField& blsP = lsP.boundaryField();
forAll(blsP, patchi)
{
const fvPatch& p = blsP[patchi].patch();
const labelUList& faceCells = p.faceCells();
// Build the d-vectors
vectorField pd(p.delta());
forAll(pd, patchFaceI)
{
dd[faceCells[patchFaceI]] +=
(1.0/magSqr(pd[patchFaceI]))*sqr(pd[patchFaceI]);
}
}
// Check for missing dimensions
// Add the missing eigenvector (such that it does not
// affect the determinant)
if (nDims == 2)
{
forAll(dd, cellI)
{
if (twoD == 0)
{
dd[cellI].xx() = 1;
}
else if (twoD == 1)
{
dd[cellI].yy() = 1;
}
else
{
dd[cellI].zz() = 1;
}
}
}
scalarField detdd(det(dd));
Info<< "max(detdd) = " << max(detdd) << nl
<< "min(detdd) = " << min(detdd) << nl
<< "average(detdd) = " << average(detdd) << endl;
label nAdaptedCells = 0;
label nAddCells = 0;
label maxNaddCells = 4*detdd.size();
additionalCellsPtr_ = new List<labelPair>(maxNaddCells);
List<labelPair>& additionalCells_ = *additionalCellsPtr_;
forAll(detdd, i)
{
label count = 0;
label oldNAddCells = nAddCells;
while (++count < 100 && detdd[i] < minDet_)
{
if (nAddCells == maxNaddCells)
{
FatalErrorIn
(
"extendedLeastSquaresVectors::"
"makeLeastSquaresVectors() const"
) << "nAddCells exceeds maxNaddCells ("
<< maxNaddCells << ")"
<< exit(FatalError);
}
labelList pointLabels = mesh.cells()[i].labels(mesh.faces());
scalar maxDetddij = 0.0;
label addCell = -1;
forAll(pointLabels, j)
{
forAll(mesh.pointCells()[pointLabels[j]], k)
{
label cellj = mesh.pointCells()[pointLabels[j]][k];
if (cellj != i)
{
vector dCij = (mesh.C()[cellj] - mesh.C()[i]);
symmTensor ddij =
dd[i] + (1.0/magSqr(dCij))*sqr(dCij);
// Add the missing eigenvector (such that it does not
// affect the determinant)
if (nDims == 2)
{
if (twoD == 0)
{
ddij.xx() = 1;
}
else if (twoD == 1)
{
ddij.yy() = 1;
}
else
{
ddij.zz() = 1;
}
}
scalar detddij = det(ddij);
if (detddij > maxDetddij)
{
addCell = cellj;
maxDetddij = detddij;
}
}
}
}
if (addCell != -1)
{
additionalCells_[nAddCells][0] = i;
additionalCells_[nAddCells++][1] = addCell;
vector dCij = mesh.C()[addCell] - mesh.C()[i];
dd[i] += (1.0/magSqr(dCij))*sqr(dCij);
// Add the missing eigenvector (such that it does not
// affect the determinant)
if (nDims == 2)
{
if (twoD == 0)
{
dd[i].xx() = 1;
}
else if (twoD == 1)
{
dd[i].yy() = 1;
}
else
{
dd[i].zz() = 1;
}
}
detdd[i] = det(dd[i]);
}
}
if (oldNAddCells < nAddCells)
{
nAdaptedCells++;
}
}
additionalCells_.setSize(nAddCells);
reduce(nAddCells, sumOp<label>());
reduce(nAdaptedCells, sumOp<label>());
if (nAddCells)
{
Info<< "max(detdd) = " << max(detdd) << nl
<< "min(detdd) = " << min(detdd) << nl
<< "average(detdd) = " << average(detdd) << nl
<< "nAdapted/nCells = "
<< scalar(nAdaptedCells)/mesh.globalData().nTotalCells() << nl
<< "nAddCells/nCells = "
<< scalar(nAddCells)/mesh.globalData().nTotalCells()
<< endl;
}
// Invert the dd tensor
const symmTensorField invDd(inv(dd));
// Revisit all faces and calculate the lsP and lsN vectors
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
vector d = C[nei] - C[own];
lsP[facei] = (1.0/magSqr(d))*(invDd[own] & d);
lsN[facei] = ((-1.0)/magSqr(d))*(invDd[nei] & d);
}
forAll(blsP, patchI)
{
fvsPatchVectorField& patchLsP = blsP[patchI];
const fvPatch& p = patchLsP.patch();
const labelUList& faceCells = p.faceCells();
// Build the d-vectors
vectorField pd(p.delta());
forAll(p, patchFaceI)
{
patchLsP[patchFaceI] =
(1.0/magSqr(pd[patchFaceI]))
*(invDd[faceCells[patchFaceI]] & pd[patchFaceI]);
}
}
additionalVectorsPtr_ = new vectorField(additionalCells_.size());
vectorField& additionalVectors_ = *additionalVectorsPtr_;
forAll(additionalCells_, i)
{
vector dCij =
mesh.C()[additionalCells_[i][1]] - mesh.C()[additionalCells_[i][0]];
additionalVectors_[i] =
(1.0/magSqr(dCij))*(invDd[additionalCells_[i][0]] & dCij);
}
if (debug)
{
Info<< "extendedLeastSquaresVectors::makeLeastSquaresVectors() :"
<< "Finished constructing least square gradient vectors"
<< endl;
}
}
const Foam::surfaceVectorField&
Foam::extendedLeastSquaresVectors::pVectors() const
{
if (!pVectorsPtr_)
{
makeLeastSquaresVectors();
}
return *pVectorsPtr_;
}
const Foam::surfaceVectorField&
Foam::extendedLeastSquaresVectors::nVectors() const
{
if (!nVectorsPtr_)
{
makeLeastSquaresVectors();
}
return *nVectorsPtr_;
}
const Foam::List<Foam::labelPair>&
Foam::extendedLeastSquaresVectors::additionalCells() const
{
if (!additionalCellsPtr_)
{
makeLeastSquaresVectors();
}
return *additionalCellsPtr_;
}
const Foam::vectorField&
Foam::extendedLeastSquaresVectors::additionalVectors() const
{
if (!additionalVectorsPtr_)
{
makeLeastSquaresVectors();
}
return *additionalVectorsPtr_;
}
bool Foam::extendedLeastSquaresVectors::movePoints()
{
deleteDemandDrivenData(pVectorsPtr_);
deleteDemandDrivenData(nVectorsPtr_);
deleteDemandDrivenData(additionalCellsPtr_);
deleteDemandDrivenData(additionalVectorsPtr_);
return true;
}
// ************************************************************************* //

190
src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -24,14 +24,13 @@ License
\*---------------------------------------------------------------------------*/
#include "leastSquaresVectors.H"
#include "surfaceFields.H"
#include "volFields.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(leastSquaresVectors, 0);
defineTypeNameAndDebug(leastSquaresVectors, 0);
}
@ -40,34 +39,7 @@ defineTypeNameAndDebug(leastSquaresVectors, 0);
Foam::leastSquaresVectors::leastSquaresVectors(const fvMesh& mesh)
:
MeshObject<fvMesh, leastSquaresVectors>(mesh),
pVectorsPtr_(NULL),
nVectorsPtr_(NULL)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::leastSquaresVectors::~leastSquaresVectors()
{
deleteDemandDrivenData(pVectorsPtr_);
deleteDemandDrivenData(nVectorsPtr_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
{
if (debug)
{
Info<< "leastSquaresVectors::makeLeastSquaresVectors() :"
<< "Constructing least square gradient vectors"
<< endl;
}
const fvMesh& mesh = mesh_;
pVectorsPtr_ = new surfaceVectorField
pVectors_
(
IOobject
(
@ -80,10 +52,8 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
),
mesh_,
dimensionedVector("zero", dimless/dimLength, vector::zero)
);
surfaceVectorField& lsP = *pVectorsPtr_;
nVectorsPtr_ = new surfaceVectorField
),
nVectors_
(
IOobject
(
@ -96,8 +66,30 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
),
mesh_,
dimensionedVector("zero", dimless/dimLength, vector::zero)
);
surfaceVectorField& lsN = *nVectorsPtr_;
)
{
calcLeastSquaresVectors();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::leastSquaresVectors::~leastSquaresVectors()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::leastSquaresVectors::calcLeastSquaresVectors()
{
if (debug)
{
Info<< "leastSquaresVectors::calcLeastSquaresVectors() :"
<< "Calculating least square gradient vectors"
<< endl;
}
const fvMesh& mesh = mesh_;
// Set local references to mesh data
const labelUList& owner = mesh_.owner();
@ -124,7 +116,8 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
}
surfaceVectorField::GeometricBoundaryField& blsP = lsP.boundaryField();
surfaceVectorField::GeometricBoundaryField& blsP =
pVectors_.boundaryField();
forAll(blsP, patchi)
{
@ -164,7 +157,7 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
const symmTensorField invDd(inv(dd));
// Revisit all faces and calculate the lsP and lsN vectors
// Revisit all faces and calculate the pVectors_ and nVectors_ vectors
forAll(owner, facei)
{
label own = owner[facei];
@ -173,8 +166,8 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
vector d = C[nei] - C[own];
scalar magSfByMagSqrd = magSf[facei]/magSqr(d);
lsP[facei] = (1 - w[facei])*magSfByMagSqrd*(invDd[own] & d);
lsN[facei] = -w[facei]*magSfByMagSqrd*(invDd[nei] & d);
pVectors_[facei] = (1 - w[facei])*magSfByMagSqrd*(invDd[own] & d);
nVectors_[facei] = -w[facei]*magSfByMagSqrd*(invDd[nei] & d);
}
forAll(blsP, patchi)
@ -214,127 +207,18 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
}
}
// For 3D meshes check the determinant of the dd tensor and switch to
// Gauss if it is less than 3
/* Currently the det(dd[celli]) criterion is incorrect: dd is weighted by Sf
if (mesh.nGeometricD() == 3)
{
label nBadCells = 0;
const cellList& cells = mesh.cells();
const scalarField& V = mesh.V();
const surfaceVectorField& Sf = mesh.Sf();
const surfaceScalarField& w = mesh.weights();
forAll(dd, celli)
{
if (det(dd[celli]) < 3)
{
nBadCells++;
const cell& c = cells[celli];
forAll(c, cellFacei)
{
label facei = c[cellFacei];
if (mesh.isInternalFace(facei))
{
scalar wf = max(min(w[facei], 0.8), 0.2);
if (celli == owner[facei])
{
lsP[facei] = (1 - wf)*Sf[facei]/V[celli];
}
else
{
lsN[facei] = -wf*Sf[facei]/V[celli];
}
}
else
{
label patchi = mesh.boundaryMesh().whichPatch(facei);
if (mesh.boundary()[patchi].size())
{
label patchFacei =
facei - mesh.boundaryMesh()[patchi].start();
if (w.boundaryField()[patchi].coupled())
{
scalar wf = max
(
min
(
w.boundaryField()[patchi][patchFacei],
0.8
),
0.2
);
lsP.boundaryField()[patchi][patchFacei] =
(1 - wf)
*Sf.boundaryField()[patchi][patchFacei]
/V[celli];
}
else
{
lsP.boundaryField()[patchi][patchFacei] =
Sf.boundaryField()[patchi][patchFacei]
/V[celli];
}
}
}
}
}
}
if (debug)
{
InfoIn("leastSquaresVectors::makeLeastSquaresVectors()")
<< "number of bad cells switched to Gauss = " << nBadCells
<< endl;
}
}
*/
if (debug)
{
Info<< "leastSquaresVectors::makeLeastSquaresVectors() :"
<< "Finished constructing least square gradient vectors"
Info<< "leastSquaresVectors::calcLeastSquaresVectors() :"
<< "Finished calculating least square gradient vectors"
<< endl;
}
}
const Foam::surfaceVectorField& Foam::leastSquaresVectors::pVectors() const
{
if (!pVectorsPtr_)
{
makeLeastSquaresVectors();
}
return *pVectorsPtr_;
}
const Foam::surfaceVectorField& Foam::leastSquaresVectors::nVectors() const
{
if (!nVectorsPtr_)
{
makeLeastSquaresVectors();
}
return *nVectorsPtr_;
}
bool Foam::leastSquaresVectors::movePoints()
{
deleteDemandDrivenData(pVectorsPtr_);
deleteDemandDrivenData(nVectorsPtr_);
calcLeastSquaresVectors();
return true;
}

22
src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -37,8 +37,7 @@ SourceFiles
#include "MeshObject.H"
#include "fvMesh.H"
#include "surfaceFieldsFwd.H"
#include "labelPair.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -56,14 +55,14 @@ class leastSquaresVectors
// Private data
//- Least-squares gradient vectors
mutable surfaceVectorField* pVectorsPtr_;
mutable surfaceVectorField* nVectorsPtr_;
surfaceVectorField pVectors_;
surfaceVectorField nVectors_;
// Private Member Functions
//- Construct Least-squares gradient vectors
void makeLeastSquaresVectors() const;
void calcLeastSquaresVectors();
public:
@ -85,11 +84,16 @@ public:
// Member functions
//- Return reference to owner least square vectors
const surfaceVectorField& pVectors() const;
const surfaceVectorField& pVectors() const
{
return pVectors_;
}
//- Return reference to neighbour least square vectors
const surfaceVectorField& nVectors() const;
const surfaceVectorField& nVectors() const
{
return nVectors_;
}
//- Delete the least square vectors when the mesh moves
virtual bool movePoints();

9
src/finiteVolume/fvMesh/extendedStencil/cellToCell/extendedCellToCellStencilTemplates.C
View File

@ -76,15 +76,14 @@ Foam::tmp
);
WeightedFieldType& wf = twf();
// cells
forAll(wf, cellI)
forAll(wf, celli)
{
const List<Type>& stField = stencilFld[cellI];
const List<WeightType>& stWeight = stencilWeights[cellI];
const List<Type>& stField = stencilFld[celli];
const List<WeightType>& stWeight = stencilWeights[celli];
forAll(stField, i)
{
wf[cellI] += stWeight[i]*stField[i];
wf[celli] += stWeight[i]*stField[i];
}
}

11
src/finiteVolume/fvMesh/fvMesh.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -38,14 +38,11 @@ License
#include "mapClouds.H"
#include "volPointInterpolation.H"
#include "extendedLeastSquaresVectors.H"
#include "extendedLeastSquaresVectors.H"
#include "leastSquaresVectors.H"
#include "CentredFitData.H"
#include "linearFitPolynomial.H"
#include "quadraticFitPolynomial.H"
#include "quadraticLinearFitPolynomial.H"
//#include "quadraticFitSnGradData.H"
#include "skewCorrectionVectors.H"
@ -132,13 +129,11 @@ void Foam::fvMesh::clearGeom()
// Things geometry dependent that are not updated.
volPointInterpolation::Delete(*this);
extendedLeastSquaresVectors::Delete(*this);
leastSquaresVectors::Delete(*this);
CentredFitData<linearFitPolynomial>::Delete(*this);
CentredFitData<quadraticFitPolynomial>::Delete(*this);
CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
skewCorrectionVectors::Delete(*this);
//quadraticFitSnGradData::Delete(*this);
// Note: should be in polyMesh::clearGeom but meshSearch not in OpenFOAM
// library
@ -154,13 +149,11 @@ void Foam::fvMesh::clearAddressing()
// Hack until proper callbacks. Below are all the fvMesh-MeshObjects.
volPointInterpolation::Delete(*this);
extendedLeastSquaresVectors::Delete(*this);
leastSquaresVectors::Delete(*this);
CentredFitData<linearFitPolynomial>::Delete(*this);
CentredFitData<quadraticFitPolynomial>::Delete(*this);
CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
skewCorrectionVectors::Delete(*this);
//quadraticFitSnGradData::Delete(*this);
centredCECCellToFaceStencilObject::Delete(*this);
centredCFCCellToFaceStencilObject::Delete(*this);
@ -718,13 +711,11 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
// Hack until proper callbacks. Below are all the fvMesh MeshObjects with a
// movePoints function.
MeshObjectMovePoints<volPointInterpolation>(*this);
MeshObjectMovePoints<extendedLeastSquaresVectors>(*this);
MeshObjectMovePoints<leastSquaresVectors>(*this);
MeshObjectMovePoints<CentredFitData<linearFitPolynomial> >(*this);
MeshObjectMovePoints<CentredFitData<quadraticFitPolynomial> >(*this);
MeshObjectMovePoints<CentredFitData<quadraticLinearFitPolynomial> >(*this);
MeshObjectMovePoints<skewCorrectionVectors>(*this);
//MeshObjectMovePoints<quadraticFitSnGradData>(*this);
return tsweptVols;
}

90
src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.C
View File

@ -24,7 +24,6 @@ License
\*---------------------------------------------------------------------------*/
#include "skewCorrectionVectors.H"
#include "surfaceFields.H"
#include "volFields.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -40,27 +39,8 @@ namespace Foam
Foam::skewCorrectionVectors::skewCorrectionVectors(const fvMesh& mesh)
:
MeshObject<fvMesh, skewCorrectionVectors>(mesh),
skew_(true),
skewCorrectionVectors_(NULL)
{}
Foam::skewCorrectionVectors::~skewCorrectionVectors()
{
deleteDemandDrivenData(skewCorrectionVectors_);
}
void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
{
if (debug)
{
Info<< "surfaceInterpolation::makeSkewCorrectionVectors() : "
<< "Constructing skew correction vectors"
<< endl;
}
skewCorrectionVectors_ = new surfaceVectorField
skew_(false),
skewCorrectionVectors_
(
IOobject
(
@ -73,8 +53,24 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
),
mesh_,
dimless
);
surfaceVectorField& SkewCorrVecs = *skewCorrectionVectors_;
)
{
calcSkewCorrectionVectors();
}
Foam::skewCorrectionVectors::~skewCorrectionVectors()
{}
void Foam::skewCorrectionVectors::calcSkewCorrectionVectors()
{
if (debug)
{
Info<< "surfaceInterpolation::calcSkewCorrectionVectors() : "
<< "Calculating skew correction vectors"
<< endl;
}
// Set local references to mesh data
const volVectorField& C = mesh_.C();
@ -92,14 +88,15 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
vector d = C[nei] - C[own];
vector Cpf = Cf[facei] - C[own];
SkewCorrVecs[facei] = Cpf - ((Sf[facei] & Cpf)/(Sf[facei] & d))*d;
skewCorrectionVectors_[facei] =
Cpf - ((Sf[facei] & Cpf)/(Sf[facei] & d))*d;
}
forAll(SkewCorrVecs.boundaryField(), patchI)
forAll(skewCorrectionVectors_.boundaryField(), patchI)
{
fvsPatchVectorField& patchSkewCorrVecs =
SkewCorrVecs.boundaryField()[patchI];
skewCorrectionVectors_.boundaryField()[patchI];
if (!patchSkewCorrVecs.coupled())
{
@ -132,21 +129,19 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
if (Sf.internalField().size())
{
skewCoeff = max(mag(SkewCorrVecs)*mesh_.deltaCoeffs()).value();
skewCoeff =
max(mag(skewCorrectionVectors_)*mesh_.deltaCoeffs()).value();
}
if (debug)
{
Info<< "surfaceInterpolation::makeSkewCorrectionVectors() : "
Info<< "surfaceInterpolation::calcSkewCorrectionVectors() : "
<< "skew coefficient = " << skewCoeff << endl;
}
//skewCoeff = 0.0;
if (skewCoeff < 1e-5)
{
skew_ = false;
deleteDemandDrivenData(skewCorrectionVectors_);
}
else
{
@ -155,43 +150,16 @@ void Foam::skewCorrectionVectors::makeSkewCorrectionVectors() const
if (debug)
{
Info<< "surfaceInterpolation::makeSkewCorrectionVectors() : "
Info<< "surfaceInterpolation::calcSkewCorrectionVectors() : "
<< "Finished constructing skew correction vectors"
<< endl;
}
}
bool Foam::skewCorrectionVectors::skew() const
{
if (skew_ == true && !skewCorrectionVectors_)
{
makeSkewCorrectionVectors();
}
return skew_;
}
const Foam::surfaceVectorField& Foam::skewCorrectionVectors::operator()() const
{
if (!skew())
{
FatalErrorIn("skewCorrectionVectors::operator()()")
<< "Cannot return skewCorrectionVectors; mesh is not skewed"
<< abort(FatalError);
}
return *skewCorrectionVectors_;
}
// Do what is neccessary if the mesh has moved
bool Foam::skewCorrectionVectors::movePoints()
{
skew_ = true;
deleteDemandDrivenData(skewCorrectionVectors_);
calcSkewCorrectionVectors();
return true;
}

22
src/finiteVolume/interpolation/surfaceInterpolation/schemes/skewCorrected/skewCorrectionVectors.H
View File

@ -37,7 +37,7 @@ SourceFiles
#include "MeshObject.H"
#include "fvMesh.H"
#include "surfaceFieldsFwd.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -57,13 +57,13 @@ class skewCorrectionVectors
// Private data
//- Is mesh skew
mutable bool skew_;
bool skew_;
//- Skew correction vectors
mutable surfaceVectorField* skewCorrectionVectors_;
surfaceVectorField skewCorrectionVectors_;
//- Construct skewness correction vectors
void makeSkewCorrectionVectors() const;
//- Calculate skewness correction vectors
void calcSkewCorrectionVectors();
public:
@ -83,12 +83,18 @@ public:
// Member functions
//- Return whether mesh is skew or not
bool skew() const;
bool skew() const
{
return skew_;
}
//- Return reference to skew vectors array
const surfaceVectorField& operator()() const;
const surfaceVectorField& operator()() const
{
return skewCorrectionVectors_;
}
//- Delete the correction vectors when the mesh moves
//- Update the correction vectors when the mesh moves
virtual bool movePoints();
};

7
src/postProcessing/functionObjects/field/fieldValues/fieldValueDelta/fieldValueDelta.C
View File

@ -39,15 +39,16 @@ namespace Foam
template<>
const char*
NamedEnum<fieldValues::fieldValueDelta::operationType, 4>::names[] =
NamedEnum<fieldValues::fieldValueDelta::operationType, 5>::names[] =
{
"add",
"subtract",
"min",
"max"
"max",
"average"
};
const NamedEnum<fieldValues::fieldValueDelta::operationType, 4>
const NamedEnum<fieldValues::fieldValueDelta::operationType, 5>
fieldValues::fieldValueDelta::operationTypeNames_;
}

4
src/regionModels/pyrolysisModels/noPyrolysis/noPyrolysis.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -49,7 +49,7 @@ void noPyrolysis::constructThermoChemistry()
{
solidChemistry_.reset
(
solidChemistryModel::New(regionMesh()).ptr()
basicSolidChemistryModel::New(regionMesh()).ptr()
);
solidThermo_.reset(&solidChemistry_->solidThermo());

6
src/regionModels/pyrolysisModels/noPyrolysis/noPyrolysis.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -37,7 +37,7 @@ SourceFiles
#include "pyrolysisModel.H"
#include "volFieldsFwd.H"
#include "solidChemistryModel.H"
#include "basicSolidChemistryModel.H"
#include "radiationModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -82,7 +82,7 @@ protected:
void constructThermoChemistry();
//- Reference to the solid chemistry model
autoPtr<solidChemistryModel> solidChemistry_;
autoPtr<basicSolidChemistryModel> solidChemistry_;
//- Reference to solid thermo
autoPtr<solidReactionThermo> solidThermo_;

6
src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -303,7 +303,7 @@ void reactingOneDim::calculateMassTransfer()
reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
:
pyrolysisModel(modelType, mesh),
solidChemistry_(solidChemistryModel::New(regionMesh())),
solidChemistry_(basicSolidChemistryModel::New(regionMesh())),
solidThermo_(solidChemistry_->solidThermo()),
radiation_(radiation::radiationModel::New(solidThermo_.T())),
rho_
@ -386,7 +386,7 @@ reactingOneDim::reactingOneDim
)
:
pyrolysisModel(modelType, mesh, dict),
solidChemistry_(solidChemistryModel::New(regionMesh())),
solidChemistry_(basicSolidChemistryModel::New(regionMesh())),
solidThermo_(solidChemistry_->solidThermo()),
radiation_(radiation::radiationModel::New(solidThermo_.T())),
rho_

6
src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -36,7 +36,7 @@ SourceFiles
#define reactingOneDim_H
#include "pyrolysisModel.H"
#include "solidChemistryModel.H"
#include "basicSolidChemistryModel.H"
#include "radiationModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -76,7 +76,7 @@ protected:
// Protected data
//- Reference to the solid chemistry model
autoPtr<solidChemistryModel> solidChemistry_;
autoPtr<basicSolidChemistryModel> solidChemistry_;
//- Reference to solid thermo
solidReactionThermo& solidThermo_;

6
src/thermophysicalModels/solidChemistryModel/Make/files
View File

@ -1,7 +1,7 @@
solidChemistryModel/solidChemistryModel.C
solidChemistryModel/solidChemistryModelNew.C
solidChemistryModel/solidChemistryModels.C
basicSolidChemistryModel/basicSolidChemistryModel.C
basicSolidChemistryModel/basicSolidChemistryModelNew.C
basicSolidChemistryModel/basicSolidChemistryModels.C
solidChemistrySolver/makeSolidChemistrySolvers.C

56
src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.C Normal file
View File

@ -0,0 +1,56 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "basicSolidChemistryModel.H"
#include "fvMesh.H"
#include "Time.H"
/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
namespace Foam
{
defineTypeNameAndDebug(basicSolidChemistryModel, 0);
defineRunTimeSelectionTable(basicSolidChemistryModel, fvMesh);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::basicSolidChemistryModel::basicSolidChemistryModel
(
const fvMesh& mesh
)
:
basicChemistryModel(mesh),
solidThermo_(solidReactionThermo::New(mesh))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::basicSolidChemistryModel::~basicSolidChemistryModel()
{}
// ************************************************************************* //

167
src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModel.H Normal file
View File

@ -0,0 +1,167 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::basicSolidChemistryModel
Description
Chemistry model for solid thermodynamics
SourceFiles
basicSolidChemistryModelI.H
basicSolidChemistryModel.C
newChemistrySolidModel.C
\*---------------------------------------------------------------------------*/
#ifndef basicSolidChemistryModel_H
#define basicSolidChemistryModel_H
#include "basicChemistryModel.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
#include "solidReactionThermo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class fvMesh;
/*---------------------------------------------------------------------------*\
class basicSolidChemistryModel Declaration
\*---------------------------------------------------------------------------*/
class basicSolidChemistryModel
:
public basicChemistryModel
{
// Private Member Functions
//- Construct as copy (not implemented)
basicSolidChemistryModel(const basicSolidChemistryModel&);
//- Disallow default bitwise assignment
void operator=(const basicSolidChemistryModel&);
protected:
// Protected data
//- Solid thermo package
autoPtr<solidReactionThermo> solidThermo_;
public:
//- Runtime type information
TypeName("basicSolidChemistryModel");
//- Declare run-time constructor selection tables
declareRunTimeSelectionTable
(
autoPtr,
basicSolidChemistryModel,
fvMesh,
(const fvMesh& mesh),
(mesh)
);
// Constructors
//- Construct from mesh
basicSolidChemistryModel(const fvMesh& mesh);
//- Selector
static autoPtr<basicSolidChemistryModel> New(const fvMesh& mesh);
//- Destructor
virtual ~basicSolidChemistryModel();
// Member Functions
//- Return access to the solid thermo package
inline solidReactionThermo& solidThermo();
//- Return const access to the solid thermo package
inline const solidReactionThermo& solidThermo() const;
//- Return total gases mass source term [kg/m3/s]
virtual tmp<DimensionedField<scalar, volMesh> > RRg() const = 0;
//- Return total solids mass source term [kg/m3/s]
virtual tmp<DimensionedField<scalar, volMesh> > RRs() const = 0;
//- Return chemical source terms for solids [kg/m3/s]
virtual const DimensionedField<scalar, volMesh>& RRs
(
const label i
) const = 0;
//- Return chemical source terms for gases [kg/m3/s]
virtual const DimensionedField<scalar, volMesh>& RRg
(
const label i
) const = 0;
//- Return sensible enthalpy for gas i [J/Kg]
virtual tmp<volScalarField> gasHs
(
const volScalarField& p,
const volScalarField& T,
const label i
) const = 0;
//- Return specie Table for gases
virtual const speciesTable& gasTable() const = 0;
//- Set reacting status of cell, cellI
virtual void setCellReacting(const label cellI, const bool active) = 0;
//- Calculates the reaction rates
virtual void calculate() = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "basicSolidChemistryModelI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

20
src/finiteVolume/finiteVolume/gradSchemes/extendedLeastSquaresGrad/extendedLeastSquaresGrads.C → src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModelI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,17 +23,19 @@ License
\*---------------------------------------------------------------------------*/
#include "fvMesh.H"
#include "extendedLeastSquaresGrad.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
inline Foam::solidReactionThermo& Foam::basicSolidChemistryModel::solidThermo()
{
return solidThermo_();
}
namespace Foam
inline const Foam::solidReactionThermo&
Foam::basicSolidChemistryModel::solidThermo() const
{
namespace fv
{
makeFvGradScheme(extendedLeastSquaresGrad)
}
return solidThermo_();
}
// ************************************************************************* //

9
src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelNew.C → src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModelNew.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,11 +23,12 @@ License
\*---------------------------------------------------------------------------*/
#include "solidChemistryModel.H"
#include "basicSolidChemistryModel.H"
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::solidChemistryModel> Foam::solidChemistryModel::New
Foam::autoPtr<Foam::basicSolidChemistryModel> Foam::basicSolidChemistryModel::
New
(
const fvMesh& mesh
)
@ -160,7 +161,7 @@ Foam::autoPtr<Foam::solidChemistryModel> Foam::solidChemistryModel::New
FatalError<< exit(FatalError);
}
return autoPtr<solidChemistryModel>(cstrIter()(mesh));
return autoPtr<basicSolidChemistryModel>(cstrIter()(mesh));
}

11
src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModels.C → src/thermophysicalModels/solidChemistryModel/basicSolidChemistryModel/basicSolidChemistryModels.C
View File

@ -22,7 +22,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
InClass
Foam::solidChemistryModel
Foam::basicSolidChemistryModel
Description
Creates solid chemistry model instances templated on the type of
@ -32,7 +32,8 @@ Description
#include "makeSolidChemistryModel.H"
#include "ODESolidChemistryModel.H"
#include "pyrolysisChemistryModel.H"
#include "basicSolidChemistryModel.H"
#include "solidChemistryModel.H"
#include "solidThermoPhysicsTypes.H"
#include "thermoPhysicsTypes.H"
@ -43,16 +44,18 @@ namespace Foam
{
makeSolidChemistryModel
(
ODESolidChemistryModel,
solidChemistryModel,
pyrolysisChemistryModel,
basicSolidChemistryModel,
hConstSolidThermoPhysics,
gasHThermoPhysics
);
makeSolidChemistryModel
(
ODESolidChemistryModel,
solidChemistryModel,
pyrolysisChemistryModel,
basicSolidChemistryModel,
hExponentialSolidThermoPhysics,
gasHThermoPhysics
);

14
src/thermophysicalModels/solidChemistryModel/makeSolidChemistryModel.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -30,6 +30,7 @@ Description
#define makeSolidChemistryModel_H
#include "addToRunTimeSelectionTable.H"
#include "solidChemistryModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -38,10 +39,19 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#define makeSolidChemistryModel(SS, Comp, SThermo, GThermo) \
#define makeSolidChemistryModel(sChemistry, SS, Comp, SThermo, GThermo) \
\
typedef SS<Comp, SThermo, GThermo> SS##Comp##SThermo##GThermo; \
\
typedef sChemistry<Comp, SThermo> sChemistryl##Comp##SThermo; \
\
defineTemplateTypeNameAndDebugWithName \
( \
sChemistryl##Comp##SThermo, \
(#sChemistry"<"#Comp"," + SThermo::typeName() + ">").c_str(), \
0 \
); \
\
defineTemplateTypeNameAndDebugWithName \
( \
SS##Comp##SThermo##GThermo, \

276
src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C → src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,73 +23,33 @@ License
\*---------------------------------------------------------------------------*/
#include "ODESolidChemistryModel.H"
#include "reactingMixture.H"
#include "pyrolysisChemistryModel.H"
#include "solidReaction.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class CompType, class SolidThermo, class GasThermo>
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
ODESolidChemistryModel
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
pyrolysisChemistryModel
(
const fvMesh& mesh
)
:
CompType(mesh),
ODE(),
Ys_(this->solidThermo().composition().Y()),
reactions_
(
dynamic_cast<const reactingMixture<SolidThermo>& >
(
this->solidThermo()
)
),
pyrolisisGases_(reactions_[0].gasSpecies()),
solidThermo_
(
dynamic_cast<const reactingMixture<SolidThermo>& >
(
this->solidThermo()
).speciesData()
),
solidChemistryModel<CompType, SolidThermo>(mesh),
pyrolisisGases_(this->reactions_[0].gasSpecies()),
gasThermo_(pyrolisisGases_.size()),
nGases_(pyrolisisGases_.size()),
nSpecie_(Ys_.size() + nGases_),
nSolids_(Ys_.size()),
nReaction_(reactions_.size()),
RRs_(nSolids_),
nSpecie_(this->Ys_.size() + nGases_),
RRg_(nGases_),
Ys0_(nSolids_),
cellCounter_(0),
reactingCells_(mesh.nCells(), true)
Ys0_(this->nSolids_),
cellCounter_(0)
{
// create the fields for the chemistry sources
forAll(RRs_, fieldI)
forAll(this->RRs_, fieldI)
{
RRs_.set
(
fieldI,
new DimensionedField<scalar, volMesh>
(
IOobject
(
"RRs." + Ys_[fieldI].name(),
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
)
);
IOobject header
(
Ys_[fieldI].name() + "0",
this->Ys_[fieldI].name() + "0",
mesh.time().timeName(),
mesh,
IOobject::NO_READ
@ -105,7 +65,7 @@ ODESolidChemistryModel
(
IOobject
(
Ys_[fieldI].name() + "0",
this->Ys_[fieldI].name() + "0",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
@ -137,7 +97,7 @@ ODESolidChemistryModel
(
IOobject
(
Ys_[fieldI].name() + "0",
this->Ys_[fieldI].name() + "0",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
@ -150,9 +110,9 @@ ODESolidChemistryModel
// Calculate inital values of Ysi0 = rho*delta*Yi
Ys0_[fieldI].internalField() =
this->solidThermo().rho()
*max(Ys_[fieldI], scalar(0.001))*mesh.V();
*max(this->Ys_[fieldI], scalar(0.001))*mesh.V();
}
}
}
forAll(RRg_, fieldI)
{
@ -190,23 +150,24 @@ ODESolidChemistryModel
);
}
Info<< "solidChemistryModel: Number of solids = " << nSolids_ << endl;
Info<< "solidChemistryModel: Number of gases = " << nGases_ << endl;
forAll(reactions_, i)
Info<< "pyrolysisChemistryModel: " << nl;
Info<< indent << "Number of solids = " << this->nSolids_ << nl;
Info<< indent << "Number of gases = " << nGases_ << nl;
forAll(this->reactions_, i)
{
Info<< indent << reactions_[i] << nl;
Info<< dynamic_cast<const solidReaction<SolidThermo>& >
(
this->reactions_[i]
) << nl;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class CompType, class SolidThermo, class GasThermo>
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
~ODESolidChemistryModel()
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
~pyrolysisChemistryModel()
{}
@ -214,7 +175,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
template<class CompType, class SolidThermo, class GasThermo>
Foam::scalarField Foam::
ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::omega
(
const scalarField& c,
const scalar T,
@ -229,9 +190,9 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
scalarField om(nEqns(), 0.0);
forAll(reactions_, i)
forAll(this->reactions_, i)
{
const Reaction<SolidThermo>& R = reactions_[i];
const Reaction<SolidThermo>& R = this->reactions_[i];
scalar omegai = omega
(
@ -242,13 +203,13 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
{
label si = R.lhs()[s].index;
om[si] -= omegai;
rhoL = solidThermo_[si].rho(p, T);
rhoL = this->solidThermo_[si].rho(p, T);
}
scalar sr = 0.0;
forAll(R.rhs(), s)
{
label si = R.rhs()[s].index;
scalar rhoR = solidThermo_[si].rho(p, T);
scalar rhoR = this->solidThermo_[si].rho(p, T);
sr = rhoR/rhoL;
om[si] += sr*omegai;
@ -260,7 +221,7 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
forAll(R.grhs(), g)
{
label gi = R.grhs()[g].index;
om[gi + nSolids_] += (1.0 - sr)*omegai;
om[gi + this->nSolids_] += (1.0 - sr)*omegai;
}
}
@ -270,7 +231,7 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
template<class CompType, class SolidThermo, class GasThermo>
Foam::scalar
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::omega
(
const Reaction<SolidThermo>& R,
const scalarField& c,
@ -314,7 +275,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
template<class CompType, class SolidThermo, class GasThermo>
void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
derivatives
(
const scalar time,
@ -331,19 +292,19 @@ derivatives
//Total mass concentration
scalar cTot = 0.0;
for (label i=0; i<nSolids_; i++)
for (label i=0; i<this->nSolids_; i++)
{
cTot += c[i];
}
scalar newCp = 0.0;
scalar newhi = 0.0;
for (label i=0; i<nSolids_; i++)
for (label i=0; i<this->nSolids_; i++)
{
scalar dYidt = dcdt[i]/cTot;
scalar Yi = c[i]/cTot;
newCp += Yi*solidThermo_[i].Cp(p, T);
newhi -= dYidt*solidThermo_[i].Hc();
newCp += Yi*this->solidThermo_[i].Cp(p, T);
newhi -= dYidt*this->solidThermo_[i].Hc();
}
scalar dTdt = newhi/newCp;
@ -356,7 +317,8 @@ derivatives
template<class CompType, class SolidThermo, class GasThermo>
void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
jacobian
(
const scalar t,
const scalarField& c,
@ -369,7 +331,7 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
scalarField c2(nSpecie_, 0.0);
for (label i=0; i<nSolids_; i++)
for (label i=0; i<this->nSolids_; i++)
{
c2[i] = max(c[i], 0.0);
}
@ -385,9 +347,9 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
// length of the first argument must be nSolids
dcdt = omega(c2, T, p);
for (label ri=0; ri<reactions_.size(); ri++)
for (label ri=0; ri<this->reactions_.size(); ri++)
{
const Reaction<SolidThermo>& R = reactions_[ri];
const Reaction<SolidThermo>& R = this->reactions_[ri];
scalar kf0 = R.kf(p, T, c2);
@ -451,96 +413,9 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
}
template<class CompType, class SolidThermo, class GasThermo>
Foam::tmp<Foam::volScalarField>
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::tc() const
{
notImplemented
(
"ODESolidChemistryModel::tc()"
);
return volScalarField::null();
}
template<class CompType, class SolidThermo, class GasThermo>
Foam::tmp<Foam::volScalarField>
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::Sh() const
{
tmp<volScalarField> tSh
(
new volScalarField
(
IOobject
(
"Sh",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
this->mesh_,
dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0),
zeroGradientFvPatchScalarField::typeName
)
);
if (this->chemistry_)
{
scalarField& Sh = tSh();
forAll(Ys_, i)
{
forAll(Sh, cellI)
{
scalar hf = solidThermo_[i].Hc();
Sh[cellI] -= hf*RRs_[i][cellI];
}
}
}
return tSh;
}
template<class CompType, class SolidThermo, class GasThermo>
Foam::tmp<Foam::volScalarField>
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::dQ() const
{
tmp<volScalarField> tdQ
(
new volScalarField
(
IOobject
(
"dQ",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh_,
dimensionedScalar("dQ", dimEnergy/dimTime, 0.0),
zeroGradientFvPatchScalarField::typeName
)
);
if (this->chemistry_)
{
volScalarField& dQ = tdQ();
dQ.dimensionedInternalField() = this->mesh_.V()*Sh()();
}
return tdQ;
}
template<class CompType, class SolidThermo, class GasThermo>
Foam::label Foam::
ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::nEqns() const
pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::nEqns() const
{
// nEqns = number of solids + gases + temperature + pressure
return (nSpecie_ + 2);
@ -548,7 +423,7 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::nEqns() const
template<class CompType, class SolidThermo, class GasThermo>
void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
calculate()
{
if (!this->chemistry_)
@ -570,10 +445,11 @@ calculate()
this->solidThermo().rho()
);
forAll(RRs_, i)
forAll(this->RRs_, i)
{
RRs_[i].field() = 0.0;
this->RRs_[i].field() = 0.0;
}
forAll(RRg_, i)
{
RRg_[i].field() = 0.0;
@ -585,28 +461,28 @@ calculate()
const scalar delta = this->mesh().V()[celli];
if (reactingCells_[celli])
if (this->reactingCells_[celli])
{
scalar rhoi = rho[celli];
scalar Ti = this->solidThermo().T()[celli];
scalar pi = this->solidThermo().p()[celli];
scalarField c(nSpecie_, 0.0);
for (label i=0; i<nSolids_; i++)
for (label i=0; i<this->nSolids_; i++)
{
c[i] = rhoi*Ys_[i][celli]*delta;
c[i] = rhoi*this->Ys_[i][celli]*delta;
}
const scalarField dcdt = omega(c, Ti, pi, true);
forAll(RRs_, i)
forAll(this->RRs_, i)
{
RRs_[i][celli] = dcdt[i]/delta;
this->RRs_[i][celli] = dcdt[i]/delta;
}
forAll(RRg_, i)
{
RRg_[i][celli] = dcdt[nSolids_ + i]/delta;
RRg_[i][celli] = dcdt[this->nSolids_ + i]/delta;
}
}
}
@ -615,7 +491,7 @@ calculate()
template<class CompType, class SolidThermo, class GasThermo>
Foam::scalar
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
(
const scalar t0,
const scalar deltaT
@ -642,9 +518,9 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
this->solidThermo().rho()
);
forAll(RRs_, i)
forAll(this->RRs_, i)
{
RRs_[i].field() = 0.0;
this->RRs_[i].field() = 0.0;
}
forAll(RRg_, i)
{
@ -654,7 +530,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
forAll(rho, celli)
{
if (reactingCells_[celli])
if (this->reactingCells_[celli])
{
cellCounter_ = celli;
@ -668,9 +544,9 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
scalar delta = this->mesh().V()[celli];
for (label i=0; i<nSolids_; i++)
for (label i=0; i<this->nSolids_; i++)
{
c[i] = rhoi*Ys_[i][celli]*delta;
c[i] = rhoi*this->Ys_[i][celli]*delta;
}
c0 = c;
@ -690,7 +566,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
scalar cTot = 0.0;
//Total mass concentration
for (label i=0; i<nSolids_; i++)
for (label i=0; i<this->nSolids_; i++)
{
cTot += c[i];
}
@ -700,13 +576,13 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
scalar invRho = 0.0;
scalarList dcdt = (c - c0)/dt;
for (label i=0; i<nSolids_; i++)
for (label i=0; i<this->nSolids_; i++)
{
scalar dYi = dcdt[i]/cTot;
scalar Yi = c[i]/cTot;
newCp += Yi*solidThermo_[i].Cp(pi, Ti);
newhi -= dYi*solidThermo_[i].Hc();
invRho += Yi/solidThermo_[i].rho(pi, Ti);
newCp += Yi*this->solidThermo_[i].Cp(pi, Ti);
newhi -= dYi*this->solidThermo_[i].Hc();
invRho += Yi/this->solidThermo_[i].rho(pi, Ti);
}
scalar dTi = (newhi/newCp)*dt;
@ -722,14 +598,14 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
deltaTMin = min(tauC, deltaTMin);
dc = c - c0;
forAll(RRs_, i)
forAll(this->RRs_, i)
{
RRs_[i][celli] = dc[i]/(deltaT*delta);
this->RRs_[i][celli] = dc[i]/(deltaT*delta);
}
forAll(RRg_, i)
{
RRg_[i][celli] = dc[nSolids_ + i]/(deltaT*delta);
RRg_[i][celli] = dc[this->nSolids_ + i]/(deltaT*delta);
}
// Update Ys0_
@ -746,7 +622,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
template<class CompType, class SolidThermo,class GasThermo>
Foam::tmp<Foam::volScalarField>
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::gasHs
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::gasHs
(
const volScalarField& p,
const volScalarField& T,
@ -785,9 +661,9 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::gasHs
}
template<class CompType, class SolidThermo,class GasThermo>
template<class CompType, class SolidThermo, class GasThermo>
Foam::scalar
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
(
scalarField &c,
const scalar T,
@ -798,7 +674,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
{
notImplemented
(
"ODESolidChemistryModel::solve"
"pyrolysisChemistryModel::solve"
"("
"scalarField&, "
"const scalar, "
@ -809,14 +685,4 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::solve
);
return (0);
}
template<class CompType, class SolidThermo,class GasThermo>
void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
setCellReacting(const label cellI, const bool active)
{
reactingCells_[cellI] = active;
}
// ************************************************************************* //

101
src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.H → src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,26 +22,24 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::ODESolidChemistryModel
Foam::pyrolysisChemistryModel
Description
Extends base chemistry model by adding a thermo package, and ODE functions.
Introduces chemistry equation system and evaluation of chemical source
terms.
Pyrolysis chemistry model. It includes gas phase in the solid
reaction.
SourceFiles
ODESolidChemistryModelI.H
ODESolidChemistryModel.C
pyrolysisChemistryModelI.H
pyrolysisChemistryModel.C
\*---------------------------------------------------------------------------*/
#ifndef ODESolidChemistryModel_H
#define ODESolidChemistryModel_H
#ifndef pyrolysisChemistryModel_H
#define pyrolysisChemistryModel_H
#include "Reaction.H"
#include "ODE.H"
#include "volFieldsFwd.H"
#include "DimensionedField.H"
#include "solidChemistryModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -52,35 +50,25 @@ namespace Foam
class fvMesh;
/*---------------------------------------------------------------------------*\
Class ODESolidChemistryModel Declaration
Class pyrolysisChemistryModel Declaration
\*---------------------------------------------------------------------------*/
template<class CompType, class SolidThermo, class GasThermo>
class ODESolidChemistryModel
class pyrolysisChemistryModel
:
public CompType,
public ODE
public solidChemistryModel<CompType, SolidThermo>
{
// Private Member Functions
//- Disallow default bitwise assignment
void operator=(const ODESolidChemistryModel&);
void operator=(const pyrolysisChemistryModel&);
protected:
//- Reference to solid mass fractions
PtrList<volScalarField>& Ys_;
//- Reactions
const PtrList<Reaction<SolidThermo> >& reactions_;
//- List of gas species present in reaction system
speciesTable pyrolisisGases_;
//- Thermodynamic data of solids
const PtrList<SolidThermo>& solidThermo_;
//- Thermodynamic data of gases
PtrList<GasThermo> gasThermo_;
@ -90,24 +78,12 @@ protected:
//- Number of components being solved by ODE
label nSpecie_;
//- Number of solid components
label nSolids_;
//- Number of solid reactions
label nReaction_;
//- List of reaction rate per solid [kg/m3/s]
PtrList<DimensionedField<scalar, volMesh> > RRs_;
//- List of reaction rate per gas [kg/m3/s]
PtrList<DimensionedField<scalar, volMesh> > RRg_;
// Protected Member Functions
//- Write access to source terms for solids
inline PtrList<DimensionedField<scalar, volMesh> >& RRs();
//- Write access to source terms for gases
inline PtrList<DimensionedField<scalar, volMesh> >& RRg();
@ -120,37 +96,25 @@ private:
//- Cell counter
label cellCounter_;
//- List of active reacting cells
List<bool> reactingCells_;
// Private members
//- Set reacting status of cell, cellI
void setCellReacting(const label cellI, const bool active);
public:
//- Runtime type information
TypeName("ODESolidChemistryModel");
TypeName("pyrolysis");
// Constructors
//- Construct from mesh
ODESolidChemistryModel(const fvMesh& mesh);
pyrolysisChemistryModel(const fvMesh& mesh);
//- Destructor
virtual ~ODESolidChemistryModel();
virtual ~pyrolysisChemistryModel();
// Member Functions
//- The reactions
inline const PtrList<Reaction<SolidThermo> >& reactions() const;
//- Thermodynamic data of gases
inline const PtrList<GasThermo>& gasThermo() const;
@ -163,9 +127,6 @@ public:
//- The number of solids
inline label nGases() const;
//- The number of reactions
inline label nReaction() const;
//- dc/dt = omega, rate of change in concentration, for each species
virtual scalarField omega
@ -198,12 +159,6 @@ public:
// Chemistry model functions
//- Return const access to the chemical source terms for solids
inline const DimensionedField<scalar, volMesh>& RRs
(
const label i
) const;
//- Return const access to the chemical source terms for gases
inline const DimensionedField<scalar, volMesh>& RRg
(
@ -213,15 +168,6 @@ public:
//- Return total gas source term
inline tmp<DimensionedField<scalar, volMesh> > RRg() const;
//- Return total solid source term
inline tmp<DimensionedField<scalar, volMesh> > RRs() const;
//- Return const access to the total source terms
inline const DimensionedField<scalar, volMesh>& RR
(
const label i
) const;
//- Return sensible enthalpy for gas i [J/Kg]
virtual tmp<volScalarField> gasHs
(
@ -232,16 +178,7 @@ public:
//- Solve the reaction system for the given start time and time
// step and return the characteristic time
virtual scalar solve(const scalar t0, const scalar deltaT);
//- Return the chemical time scale
virtual tmp<volScalarField> tc() const;
//- Return source for enthalpy equation [kg/m/s3]
virtual tmp<volScalarField> Sh() const;
//- Return the heat release, i.e. enthalpy/sec [m2/s3]
virtual tmp<volScalarField> dQ() const;
virtual scalar solve(const scalar t0, const scalar deltaT) ;
// ODE functions (overriding abstract functions in ODE.H)
@ -281,12 +218,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
# include "ODESolidChemistryModelI.H"
# include "pyrolysisChemistryModelI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "ODESolidChemistryModel.C"
# include "pyrolysisChemistryModel.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

102
src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModelI.H → src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModelI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -24,37 +24,21 @@ License
\*---------------------------------------------------------------------------*/
#include "volFields.H"
#include "zeroGradientFvPatchFields.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CompType, class SolidThermo, class GasThermo>
inline Foam::PtrList<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >&
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRs()
{
return RRs_;
}
template<class CompType, class SolidThermo, class GasThermo>
inline Foam::PtrList<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >&
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg()
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::RRg()
{
return RRg_;
}
template<class CompType, class SolidThermo, class GasThermo>
inline const Foam::PtrList<Foam::Reaction<SolidThermo> >&
Foam::ODESolidChemistryModel<CompType, SolidThermo,GasThermo>::reactions() const
{
return reactions_;
}
template<class CompType, class SolidThermo, class GasThermo>
inline const Foam::PtrList<GasThermo>&
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
gasThermo() const
{
return gasThermo_;
@ -63,7 +47,8 @@ gasThermo() const
template<class CompType, class SolidThermo, class GasThermo>
inline const Foam::speciesTable&
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::gasTable() const
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
gasTable() const
{
return pyrolisisGases_;
}
@ -71,35 +56,16 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::gasTable() const
template<class CompType, class SolidThermo, class GasThermo>
inline Foam::label
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::nSpecie() const
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
nSpecie() const
{
return nSpecie_;
}
template<class CompType, class SolidThermo, class GasThermo>
inline Foam::label
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
nReaction() const
{
return nReaction_;
}
template<class CompType, class SolidThermo, class GasThermo>
inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRs
(
const label i
) const
{
return RRs_[i];
}
template<class CompType, class SolidThermo, class GasThermo>
inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::RRg
(
const label i
) const
@ -110,7 +76,8 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg
template<class CompType, class SolidThermo, class GasThermo>
inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg() const
Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
RRg() const
{
tmp<DimensionedField<scalar, volMesh> > tRRg
(
@ -139,51 +106,4 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg() const
}
return tRRg;
}
template<class CompType, class SolidThermo, class GasThermo>
inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRs() const
{
tmp<DimensionedField<scalar, volMesh> > tRRs
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
"RRs",
this->time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
this->mesh(),
dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
)
);
if (this->chemistry_)
{
DimensionedField<scalar, volMesh>& RRs = tRRs();
for (label i=0; i < nSolids_; i++)
{
RRs += RRs_[i];
}
}
return tRRs;
}
template<class CompType, class SolidThermo, class GasThermo>
inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RR
(
const label i
) const
{
notImplemented("ODESolidChemistryModel::RR(const label)");
return (DimensionedField<scalar, volMesh>::null());
}
// ************************************************************************* //

192
src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -24,33 +24,193 @@ License
\*---------------------------------------------------------------------------*/
#include "solidChemistryModel.H"
#include "fvMesh.H"
#include "Time.H"
/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
namespace Foam
{
defineTypeNameAndDebug(solidChemistryModel, 0);
defineRunTimeSelectionTable(solidChemistryModel, fvMesh);
}
#include "reactingMixture.H"
#include "zeroGradientFvPatchFields.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::solidChemistryModel::solidChemistryModel
template<class CompType, class SolidThermo>
Foam::solidChemistryModel<CompType, SolidThermo>::
solidChemistryModel
(
const fvMesh& mesh
)
:
basicChemistryModel(mesh),
solidThermo_(solidReactionThermo::New(mesh))
{}
CompType(mesh),
ODE(),
Ys_(this->solidThermo().composition().Y()),
reactions_
(
dynamic_cast<const reactingMixture<SolidThermo>& >
(
this->solidThermo()
)
),
solidThermo_
(
dynamic_cast<const reactingMixture<SolidThermo>& >
(
this->solidThermo()
).speciesData()
),
nSolids_(Ys_.size()),
nReaction_(reactions_.size()),
RRs_(nSolids_),
reactingCells_(mesh.nCells(), true)
{
// create the fields for the chemistry sources
forAll(RRs_, fieldI)
{
RRs_.set
(
fieldI,
new DimensionedField<scalar, volMesh>
(
IOobject
(
"RRs." + Ys_[fieldI].name(),
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
)
);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::solidChemistryModel::~solidChemistryModel()
template<class CompType, class SolidThermo>
Foam::solidChemistryModel<CompType, SolidThermo>::
~solidChemistryModel()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CompType, class SolidThermo>
Foam::tmp<Foam::volScalarField>
Foam::solidChemistryModel<CompType, SolidThermo>::tc() const
{
notImplemented
(
"solidChemistryModel::tc()"
);
return volScalarField::null();
}
template<class CompType, class SolidThermo>
Foam::tmp<Foam::volScalarField>
Foam::solidChemistryModel<CompType, SolidThermo>::Sh() const
{
tmp<volScalarField> tSh
(
new volScalarField
(
IOobject
(
"Sh",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
this->mesh_,
dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0),
zeroGradientFvPatchScalarField::typeName
)
);
if (this->chemistry_)
{
scalarField& Sh = tSh();
forAll(Ys_, i)
{
forAll(Sh, cellI)
{
scalar hf = solidThermo_[i].Hc();
Sh[cellI] -= hf*RRs_[i][cellI];
}
}
}
return tSh;
}
template<class CompType, class SolidThermo>
Foam::tmp<Foam::volScalarField>
Foam::solidChemistryModel<CompType, SolidThermo>::dQ() const
{
tmp<volScalarField> tdQ
(
new volScalarField
(
IOobject
(
"dQ",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh_,
dimensionedScalar("dQ", dimEnergy/dimTime, 0.0),
zeroGradientFvPatchScalarField::typeName
)
);
if (this->chemistry_)
{
volScalarField& dQ = tdQ();
dQ.dimensionedInternalField() = this->mesh_.V()*Sh()();
}
return tdQ;
}
template<class CompType, class SolidThermo>
Foam::scalar Foam::solidChemistryModel<CompType, SolidThermo>::solve
(
scalarField &c,
const scalar T,
const scalar p,
const scalar t0,
const scalar dt
) const
{
notImplemented
(
"solidChemistryModel::solve"
"("
"scalarField&, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar"
") const"
);
return (0);
}
template<class CompType, class SolidThermo>
void Foam::solidChemistryModel<CompType, SolidThermo>::setCellReacting
(
const label cellI,
const bool active
)
{
reactingCells_[cellI] = active;
}
// ************************************************************************* //

195
src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,23 +25,24 @@ Class
Foam::solidChemistryModel
Description
Chemistry model for solid thermodynamics
Extends base solid chemistry model by adding a thermo package, and ODE
functions.
Introduces chemistry equation system and evaluation of chemical source
terms.
SourceFiles
solidChemistryModelI.H
solidChemistryModel.C
newChemistrySolidModel.C
\*---------------------------------------------------------------------------*/
#ifndef solidChemistryModel_H
#define solidChemistryModel_H
#include "basicChemistryModel.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
#include "solidReactionThermo.H"
#include "Reaction.H"
#include "ODE.H"
#include "volFieldsFwd.H"
#include "DimensionedField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -52,45 +53,57 @@ namespace Foam
class fvMesh;
/*---------------------------------------------------------------------------*\
class solidChemistryModel Declaration
Class solidChemistryModel Declaration
\*---------------------------------------------------------------------------*/
template<class CompType, class SolidThermo>
class solidChemistryModel
:
public basicChemistryModel
public CompType,
public ODE
{
// Private Member Functions
//- Construct as copy (not implemented)
solidChemistryModel(const solidChemistryModel&);
//- Disallow default bitwise assignment
void operator=(const solidChemistryModel&);
protected:
// Protected data
//- Reference to solid mass fractions
PtrList<volScalarField>& Ys_;
//- Solid thermo package
autoPtr<solidReactionThermo> solidThermo_;
//- Reactions
const PtrList<Reaction<SolidThermo> >& reactions_;
//- Thermodynamic data of solids
const PtrList<SolidThermo>& solidThermo_;
//- Number of solid components
label nSolids_;
//- Number of solid reactions
label nReaction_;
//- List of reaction rate per solid [kg/m3/s]
PtrList<DimensionedField<scalar, volMesh> > RRs_;
//- List of active reacting cells
List<bool> reactingCells_;
// Protected Member Functions
//- Write access to source terms for solids
inline PtrList<DimensionedField<scalar, volMesh> >& RRs();
//- Set reacting status of cell, cellI
void setCellReacting(const label cellI, const bool active);
public:
//- Runtime type information
TypeName("solid");
//- Declare run-time constructor selection tables
declareRunTimeSelectionTable
(
autoPtr,
solidChemistryModel,
fvMesh,
(const fvMesh& mesh),
(mesh)
);
TypeName("solidChemistryModel");
// Constructors
@ -99,56 +112,108 @@ public:
solidChemistryModel(const fvMesh& mesh);
//- Selector
static autoPtr<solidChemistryModel> New(const fvMesh& mesh);
//- Destructor
virtual ~solidChemistryModel();
// Member Functions
//- Return access to the solid thermo package
inline solidReactionThermo& solidThermo();
//- The reactions
inline const PtrList<Reaction<SolidThermo> >& reactions() const;
//- Return const access to the solid thermo package
inline const solidReactionThermo& solidThermo() const;
//- The number of reactions
inline label nReaction() const;
//- Return total gases mass source term [kg/m3/s]
virtual tmp<DimensionedField<scalar, volMesh> > RRg() const = 0;
//- Return total solids mass source term [kg/m3/s]
virtual tmp<DimensionedField<scalar, volMesh> > RRs() const = 0;
//- Return chemical source terms for solids [kg/m3/s]
virtual const DimensionedField<scalar, volMesh>& RRs
//- dc/dt = omega, rate of change in concentration, for each species
virtual scalarField omega
(
const label i
const scalarField& c,
const scalar T,
const scalar p,
const bool updateC0 = false
) const = 0;
//- Return chemical source terms for gases [kg/m3/s]
virtual const DimensionedField<scalar, volMesh>& RRg
//- Return the reaction rate for reaction r and the reference
// species and charateristic times
virtual scalar omega
(
const label i
const Reaction<SolidThermo>& r,
const scalarField& c,
const scalar T,
const scalar p,
scalar& pf,
scalar& cf,
label& lRef,
scalar& pr,
scalar& cr,
label& rRef
) const = 0;
//- Return sensible enthalpy for gas i [J/Kg]
virtual tmp<volScalarField> gasHs
(
const volScalarField& p,
const volScalarField& T,
const label i
) const = 0;
//- Return specie Table for gases
virtual const speciesTable& gasTable() const = 0;
//- Set reacting status of cell, cellI
virtual void setCellReacting(const label cellI, const bool active) = 0;
//- Calculates the reaction rates
virtual void calculate() = 0;
// Chemistry model functions
//- Return const access to the chemical source terms for solids
inline const DimensionedField<scalar, volMesh>& RRs
(
const label i
) const;
//- Return total solid source term
inline tmp<DimensionedField<scalar, volMesh> > RRs() const;
//- Return const access to the total source terms
inline const DimensionedField<scalar, volMesh>& RR
(
const label i
) const;
//- Solve the reaction system for the given start time and time
// step and return the characteristic time
virtual scalar solve(const scalar t0, const scalar deltaT) = 0;
//- Return the chemical time scale
virtual tmp<volScalarField> tc() const;
//- Return source for enthalpy equation [kg/m/s3]
virtual tmp<volScalarField> Sh() const;
//- Return the heat release, i.e. enthalpy/sec [m2/s3]
virtual tmp<volScalarField> dQ() const;
// ODE functions (overriding abstract functions in ODE.H)
//- Number of ODE's to solve
virtual label nEqns() const = 0;
virtual void derivatives
(
const scalar t,
const scalarField& c,
scalarField& dcdt
) const = 0;
virtual void jacobian
(
const scalar t,
const scalarField& c,
scalarField& dcdt,
scalarSquareMatrix& dfdc
) const = 0;
virtual scalar solve
(
scalarField &c,
const scalar T,
const scalar p,
const scalar t0,
const scalar dt
) const = 0;
};
@ -158,7 +223,13 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "solidChemistryModelI.H"
# include "solidChemistryModelI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "solidChemistryModel.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

81
src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModelI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,18 +23,87 @@ License
\*---------------------------------------------------------------------------*/
#include "volFields.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline Foam::solidReactionThermo& Foam::solidChemistryModel::solidThermo()
template<class CompType, class SolidThermo>
inline Foam::PtrList<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >&
Foam::solidChemistryModel<CompType, SolidThermo>::RRs()
{
return solidThermo_();
return RRs_;
}
template<class CompType, class SolidThermo>
inline const Foam::PtrList<Foam::Reaction<SolidThermo> >&
Foam::solidChemistryModel<CompType, SolidThermo>::reactions() const
{
return reactions_;
}
inline const Foam::solidReactionThermo&
Foam::solidChemistryModel::solidThermo() const
template<class CompType, class SolidThermo>
inline Foam::label
Foam::solidChemistryModel<CompType, SolidThermo>::
nReaction() const
{
return solidThermo_();
return nReaction_;
}
template<class CompType, class SolidThermo>
inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::solidChemistryModel<CompType, SolidThermo>::RRs
(
const label i
) const
{
return RRs_[i];
}
template<class CompType, class SolidThermo>
inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::solidChemistryModel<CompType, SolidThermo>::RRs() const
{
tmp<DimensionedField<scalar, volMesh> > tRRs
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
"RRs",
this->time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
this->mesh(),
dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
)
);
if (this->chemistry_)
{
DimensionedField<scalar, volMesh>& RRs = tRRs();
for (label i=0; i < nSolids_; i++)
{
RRs += RRs_[i];
}
}
return tRRs;
}
template<class CompType, class SolidThermo>
inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::solidChemistryModel<CompType, SolidThermo>::RR
(
const label i
) const
{
notImplemented("solidChemistryModel::RR(const label)");
return (DimensionedField<scalar, volMesh>::null());
}

14
src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolverType.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -39,15 +39,15 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#define makeSolidChemistrySolverType(SS, Comp, SThermo, GThermo) \
#define makeSolidChemistrySolverType(SS, Schem, Comp, SThermo, GThermo) \
\
typedef SS<ODESolidChemistryModel<Comp, SThermo, GThermo> > \
SS##Comp##SThermo##GThermo; \
typedef SS<Schem<Comp, SThermo, GThermo> > \
SS##Schem##Comp##SThermo##GThermo; \
\
defineTemplateTypeNameAndDebugWithName \
( \
SS##Comp##SThermo##GThermo, \
(#SS"<" + word(Comp::typeName_()) \
SS##Schem##Comp##SThermo##GThermo, \
(#SS"<" + word(Schem::typeName_()) \
+ "," + SThermo::typeName() + "," + GThermo::typeName() + ">").c_str(), \
0 \
); \
@ -55,7 +55,7 @@ namespace Foam
addToRunTimeSelectionTable \
( \
Comp, \
SS##Comp##SThermo##GThermo, \
SS##Schem##Comp##SThermo##GThermo, \
fvMesh \
);

10
src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolvers.C
View File

@ -27,8 +27,8 @@ License
#include "solidThermoPhysicsTypes.H"
#include "thermoPhysicsTypes.H"
#include "ODESolidChemistryModel.H"
#include "solidChemistryModel.H"
#include "pyrolysisChemistryModel.H"
#include "basicSolidChemistryModel.H"
#include "ode.H"
@ -39,7 +39,8 @@ namespace Foam
makeSolidChemistrySolverType
(
ode,
solidChemistryModel,
pyrolysisChemistryModel,
basicSolidChemistryModel,
hConstSolidThermoPhysics,
gasHThermoPhysics
)
@ -47,7 +48,8 @@ namespace Foam
makeSolidChemistrySolverType
(
ode,
solidChemistryModel,
pyrolysisChemistryModel,
basicSolidChemistryModel,
hExponentialSolidThermoPhysics,
gasHThermoPhysics
)

20
src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2013-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -150,12 +150,18 @@ Foam::string Foam::solidReaction<ReactionThermo>::solidReactionStr
) const
{
this->reactionStrLeft(reaction);
reaction << " + ";
solidReactionStrLeft(reaction);
if (glhs().size() > 0)
{
reaction << " + ";
solidReactionStrLeft(reaction);
}
reaction << " = ";
this->reactionStrRight(reaction);
reaction << " + ";
solidReactionStrRight(reaction);
if (grhs().size() > 0)
{
reaction << " + ";
solidReactionStrRight(reaction);
}
return reaction.str();
}
@ -169,8 +175,6 @@ void Foam::solidReaction<ReactionThermo>::solidReactionStrLeft
{
for (label i = 0; i < glhs().size(); ++i)
{
reaction << " + ";
if (i > 0)
{
reaction << " + ";
@ -197,8 +201,6 @@ void Foam::solidReaction<ReactionThermo>::solidReactionStrRight
for (label i = 0; i < grhs().size(); ++i)
{
reaction << " + ";
if (i > 0)
{
reaction << " + ";