Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
Henry
2012-08-02 15:30:17 +01:00
64 changed files with 2445 additions and 769 deletions

View File

@ -97,6 +97,16 @@ void Foam::polyMesh::clearAddressing()
geometricD_ = Vector<label>::zero;
solutionD_ = Vector<label>::zero;
// Update zones
pointZones_.clearAddressing();
faceZones_.clearAddressing();
cellZones_.clearAddressing();
// Remove the stored tet base points
tetBasePtIsPtr_.clear();
// Remove the cell tree
cellTreePtr_.clear();
pointMesh::Delete(*this);
}

View File

@ -0,0 +1,52 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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/>.
\*---------------------------------------------------------------------------*/
#ifndef DataEntryFws_H
#define DataEntryFws_H
#include "DataEntry.H"
#include "vector.H"
#include "symmTensor.H"
#include "sphericalTensor.H"
#include "tensor.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
typedef DataEntry<label> labelDataEntry;
typedef DataEntry<scalar> scalarDataEntry;
typedef DataEntry<vector> vectorDataEntry;
typedef DataEntry<symmTensor> symmTensorDataEntry;
typedef DataEntry<sphericalTensor> sphericalTensorDataEntry;
typedef DataEntry<tensor> tensorDataEntry;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -25,14 +25,14 @@ License
#include "polynomial.H"
#include "Time.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(polynomial, 0);
DataEntry<scalar>::adddictionaryConstructorToTable<polynomial>
addpolynomialConstructorToTable_;
addToRunTimeSelectionTable(scalarDataEntry, polynomial, dictionary);
}
@ -40,7 +40,7 @@ namespace Foam
Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
:
DataEntry<scalar>(entryName),
scalarDataEntry(entryName),
coeffs_(),
canIntegrate_(true),
dimensions_(dimless)
@ -52,15 +52,17 @@ Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
is.putBack(firstToken);
if (firstToken == token::BEGIN_SQR)
{
is >> this->dimensions_;
is >> this->dimensions_;
}
is >> coeffs_;
if (!coeffs_.size())
{
FatalErrorIn("Foam::polynomial::polynomial(const word&, dictionary&)")
<< "polynomial coefficients for entry " << this->name_
FatalErrorIn
(
"Foam::polynomial::polynomial(const word&, const dictionary&)"
) << "polynomial coefficients for entry " << this->name_
<< " are invalid (empty)" << nl << exit(FatalError);
}
@ -77,8 +79,10 @@ Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
{
if (!canIntegrate_)
{
WarningIn("Foam::polynomial::polynomial(const word&, dictionary&)")
<< "Polynomial " << this->name_ << " cannot be integrated"
WarningIn
(
"Foam::polynomial::polynomial(const word&, const dictionary&)"
) << "Polynomial " << this->name_ << " cannot be integrated"
<< endl;
}
}
@ -91,7 +95,7 @@ Foam::polynomial::polynomial
const List<Tuple2<scalar, scalar> >& coeffs
)
:
DataEntry<scalar>(entryName),
scalarDataEntry(entryName),
coeffs_(coeffs),
canIntegrate_(true),
dimensions_(dimless)
@ -101,7 +105,7 @@ Foam::polynomial::polynomial
FatalErrorIn
(
"Foam::polynomial::polynomial"
"(const word&, const List<Tuple2<scalar, scalar> >&&)"
"(const word&, const List<Tuple2<scalar, scalar> >&)"
) << "polynomial coefficients for entry " << this->name_
<< " are invalid (empty)" << nl << exit(FatalError);
}
@ -122,7 +126,7 @@ Foam::polynomial::polynomial
WarningIn
(
"Foam::polynomial::polynomial"
"(const word&, const List<Tuple2<scalar, scalar> >&&)"
"(const word&, const List<Tuple2<scalar, scalar> >&)"
) << "Polynomial " << this->name_ << " cannot be integrated"
<< endl;
}
@ -132,7 +136,7 @@ Foam::polynomial::polynomial
Foam::polynomial::polynomial(const polynomial& poly)
:
DataEntry<scalar>(poly),
scalarDataEntry(poly),
coeffs_(poly.coeffs_),
canIntegrate_(poly.canIntegrate_),
dimensions_(poly.dimensions_)
@ -201,7 +205,8 @@ Foam::dimensioned<Foam::scalar> Foam::polynomial::dimValue
Foam::dimensioned<Foam::scalar> Foam::polynomial::dimIntegrate
(
const scalar x1, const scalar x2
const scalar x1,
const scalar x2
) const
{
return dimensioned<scalar>
@ -212,4 +217,5 @@ Foam::dimensioned<Foam::scalar> Foam::polynomial::dimIntegrate
);
}
// ************************************************************************* //

View File

@ -48,6 +48,7 @@ SourceFiles
#include "DataEntry.H"
#include "Tuple2.H"
#include "dimensionSet.H"
#include "DataEntryFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -70,7 +71,7 @@ Ostream& operator<<
class polynomial
:
public DataEntry<scalar>
public scalarDataEntry
{
// Private data
@ -107,9 +108,9 @@ public:
polynomial(const polynomial& poly);
//- Construct and return a clone
virtual tmp<DataEntry<scalar> > clone() const
virtual tmp<scalarDataEntry> clone() const
{
return tmp<DataEntry<scalar> >(new polynomial(*this));
return tmp<scalarDataEntry>(new polynomial(*this));
}

View File

@ -10,6 +10,7 @@ rhoCombustionModel/rhoCombustionModel/rhoCombustionModelNew.C
rhoCombustionModel/rhoThermoCombustion/rhoThermoCombustion.C
rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C
diffusion/diffusions.C
infinitelyFastChemistry/infinitelyFastChemistrys.C
PaSR/PaSRs.C

View File

@ -138,7 +138,7 @@ void Foam::combustionModels::PaSR<Type>::correct()
template<class Type>
Foam::tmp<Foam::fvScalarMatrix>
Foam::combustionModels::PaSR<Type>::R(const volScalarField& Y) const
Foam::combustionModels::PaSR<Type>::R(volScalarField& Y) const
{
tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));

View File

@ -103,7 +103,7 @@ public:
virtual void correct();
//- Fuel consumption rate matrix.
virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
//- Heat release rate calculated from fuel consumption rate matrix
virtual tmp<volScalarField> dQ() const;

View File

@ -135,7 +135,7 @@ public:
virtual void correct() = 0;
//- Fuel consumption rate matrix, i.e. source term for fuel equation
virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const = 0;
virtual tmp<fvScalarMatrix> R(volScalarField& Y) const = 0;
//- Heat release rate calculated from fuel consumption rate matrix
virtual tmp<volScalarField> dQ() const = 0;

View File

@ -0,0 +1,106 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
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 "diffusion.H"
#include "fvcGrad.H"
namespace Foam
{
namespace combustionModels
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class CombThermoType, class ThermoType>
diffusion<CombThermoType, ThermoType>::diffusion
(
const word& modelType, const fvMesh& mesh
)
:
singleStepCombustion<CombThermoType, ThermoType>(modelType, mesh),
C_(readScalar(this->coeffs().lookup("C"))),
oxidantName_(this->coeffs().template lookupOrDefault<word>("oxidant", "O2"))
{}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
template<class CombThermoType, class ThermoType>
diffusion<CombThermoType, ThermoType>::~diffusion()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class CombThermoType, class ThermoType>
void diffusion<CombThermoType, ThermoType>::correct()
{
this->wFuel_ ==
dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0);
if (this->active())
{
this->singleMixturePtr_->fresCorrect();
const label fuelI = this->singleMixturePtr_->fuelIndex();
const volScalarField& YFuel =
this->thermoPtr_->composition().Y()[fuelI];
if (this->thermoPtr_->composition().contains(oxidantName_))
{
const volScalarField& YO2 =
this->thermoPtr_->composition().Y(oxidantName_);
this->wFuel_ ==
C_*this->turbulence().muEff()
*mag(fvc::grad(YFuel) & fvc::grad(YO2))
*pos(YFuel)*pos(YO2);
}
}
}
template<class CombThermoType, class ThermoType>
bool diffusion<CombThermoType, ThermoType>::read()
{
if (singleStepCombustion<CombThermoType, ThermoType>::read())
{
this->coeffs().lookup("C") >> C_ ;
this->coeffs().readIfPresent("oxidant", oxidantName_);
return true;
}
else
{
return false;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace combustionModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -0,0 +1,122 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
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::combustionModels::diffusion
Description
Simple diffusion-based combustion model based on the principle mixed is
burnt. Additional parameter C is used to distribute the heat release rate
in time.
SourceFiles
diffusion.C
\*---------------------------------------------------------------------------*/
#ifndef diffusion_H
#define diffusion_H
#include "singleStepCombustion.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace combustionModels
{
/*---------------------------------------------------------------------------*\
Class diffusion Declaration
\*---------------------------------------------------------------------------*/
template<class CombThermoType, class ThermoType>
class diffusion
:
public singleStepCombustion<CombThermoType, ThermoType>
{
// Private data
//- Model constant
scalar C_;
//- Name of oxidant - default is "O2"
word oxidantName_;
// Private Member Functions
//- Disallow copy construct
diffusion(const diffusion&);
//- Disallow default bitwise assignment
void operator=(const diffusion&);
public:
//- Runtime type information
TypeName("diffusion");
// Constructors
//- Construct from components
diffusion(const word& modelType, const fvMesh& mesh);
//- Destructor
virtual ~diffusion();
// Member Functions
// Evolution
//- Correct combustion rate
virtual void correct();
// I-O
//- Update properties
virtual bool read();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace combustionModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "diffusion.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,74 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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 "makeCombustionTypes.H"
#include "thermoPhysicsTypes.H"
#include "psiThermoCombustion.H"
#include "rhoThermoCombustion.H"
#include "diffusion.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace combustionModels
{
makeCombustionTypesThermo
(
diffusion,
psiThermoCombustion,
gasThermoPhysics,
psiCombustionModel
);
makeCombustionTypesThermo
(
diffusion,
psiThermoCombustion,
constGasThermoPhysics,
psiCombustionModel
);
makeCombustionTypesThermo
(
diffusion,
rhoThermoCombustion,
gasThermoPhysics,
rhoCombustionModel
);
makeCombustionTypesThermo
(
diffusion,
rhoThermoCombustion,
constGasThermoPhysics,
rhoCombustionModel
);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -59,7 +59,7 @@ template<class CombThermoType>
Foam::tmp<Foam::fvScalarMatrix>
Foam::combustionModels::noCombustion<CombThermoType>::R
(
const volScalarField& Y
volScalarField& Y
) const
{
tmp<fvScalarMatrix> tSu

View File

@ -84,8 +84,8 @@ public:
//- Correct combustion rate
virtual void correct();
//- Fuel consumption rate matrix.
virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
//- Fuel consumption rate matrix
virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
//- Heat release rate calculated from fuel consumption rate matrix
virtual tmp<volScalarField> dQ() const;

View File

@ -53,7 +53,8 @@ singleStepCombustion<CombThermoType, ThermoType>::singleStepCombustion
),
this->mesh(),
dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
)
),
semiImplicit_(readBool(this->coeffs_.lookup("semiImplicit")))
{
if (isA<singleStepReactingMixture<ThermoType> >(this->thermo()))
{
@ -79,6 +80,15 @@ singleStepCombustion<CombThermoType, ThermoType>::singleStepCombustion
<< "Please select a thermo package based on "
<< "singleStepReactingMixture" << exit(FatalError);
}
if (semiImplicit_)
{
Info<< "Combustion mode: semi-implicit" << endl;
}
else
{
Info<< "Combustion mode: explicit" << endl;
}
}
@ -95,17 +105,28 @@ template<class CombThermoType, class ThermoType>
Foam::tmp<Foam::fvScalarMatrix>
singleStepCombustion<CombThermoType, ThermoType>::R
(
const volScalarField& Y
volScalarField& Y
) const
{
const label specieI = this->thermoPtr_->composition().species()[Y.name()];
const volScalarField wSpecie
volScalarField wSpecie
(
wFuel_*singleMixturePtr_->specieStoichCoeffs()[specieI]
);
return wSpecie + fvm::Sp(0.0*wSpecie, Y);
if (semiImplicit_)
{
const label fNorm = singleMixturePtr_->specieProd()[specieI];
const volScalarField fres = singleMixturePtr_->fres(specieI);
wSpecie /= max(fNorm*(Y - fres), 1e-2);
return -fNorm*wSpecie*fres + fNorm*fvm::Sp(wSpecie, Y);
}
else
{
return wSpecie + fvm::Sp(0.0*wSpecie, Y);
}
}
@ -114,7 +135,8 @@ Foam::tmp<Foam::volScalarField>
singleStepCombustion< CombThermoType, ThermoType>::Sh() const
{
const label fuelI = singleMixturePtr_->fuelIndex();
const volScalarField& YFuel = this->thermoPtr_->composition().Y(fuelI);
volScalarField& YFuel =
const_cast<volScalarField&>(this->thermoPtr_->composition().Y(fuelI));
return -singleMixturePtr_->qFuel()*(R(YFuel) & YFuel);
}

View File

@ -71,6 +71,9 @@ protected:
//- Fuel consumption rate
volScalarField wFuel_;
//- Semi-implicit (true) or explicit (false) treatment
bool semiImplicit_;
public:
@ -89,7 +92,7 @@ public:
// Evolution
//- Fuel consumption rate matrix
virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
//- Heat release rate calculated from fuel consumption rate matrix
virtual tmp<volScalarField> dQ() const;

View File

@ -1,4 +1,5 @@
/* Coal parcel and sub-models */
coalParcel/makeCoalParcelSubmodels.C
coalCloudList/coalCloudList.C
LIB = $(FOAM_LIBBIN)/libcoalCombustion

View File

@ -0,0 +1,92 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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 "coalCloudList.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::coalCloudList::coalCloudList
(
const volScalarField& rho,
const volVectorField& U,
const dimensionedVector& g,
const SLGThermo& slgThermo
)
:
PtrList<coalCloud>(),
mesh_(rho.mesh())
{
IOdictionary props
(
IOobject
(
"coalCloudList",
mesh_.time().constant(),
mesh_,
IOobject::MUST_READ
)
);
const wordHashSet cloudNames(wordList(props.lookup("clouds")));
setSize(cloudNames.size());
label i = 0;
forAllConstIter(wordHashSet, cloudNames, iter)
{
const word& name = iter.key();
Info<< "creating cloud: " << name << endl;
set
(
i++,
new coalCloud
(
name,
rho,
U,
g,
slgThermo
)
);
Info<< endl;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::coalCloudList::evolve()
{
forAll(*this, i)
{
operator[](i).evolve();
}
}
// ************************************************************************* //

View File

@ -0,0 +1,126 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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/>.
\*---------------------------------------------------------------------------*/
#ifndef coalCloudList_H
#define coalCloudList_H
#include "coalCloud.H"
#include "volFieldsFwd.H"
#include "fvMatricesFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class coalCloudList Declaration
\*---------------------------------------------------------------------------*/
class coalCloudList
:
public PtrList<coalCloud>
{
private:
//- Reference to the mesh
const fvMesh& mesh_;
public:
// Constructor
coalCloudList
(
const volScalarField& rho,
const volVectorField& U,
const dimensionedVector& g,
const SLGThermo& slgThermo
);
// Member Functions
// Evolution
//- Evolve the cloud collection
void evolve();
// Source terms
//- Return const reference to momentum source
inline tmp<DimensionedField<vector, volMesh> > UTrans() const;
//- Return tmp momentum source term
inline tmp<fvVectorMatrix> SU(volVectorField& U) const;
//- Sensible enthalpy transfer [J/kg]
inline tmp<DimensionedField<scalar, volMesh> > hsTrans() const;
//- Return sensible enthalpy source term [J/kg/m3/s]
inline tmp<fvScalarMatrix> Sh(volScalarField& hs) const;
//- Return mass source term for specie i - specie eqn
inline tmp<fvScalarMatrix> SYi
(
const label i,
volScalarField& Yi
) const;
//- Return total mass transfer [kg/m3]
inline tmp<DimensionedField<scalar, volMesh> > rhoTrans() const;
//- Return tmp total mass source for carrier phase
// - fully explicit
inline tmp<DimensionedField<scalar, volMesh> > Srho() const;
//- Return tmp total mass source for carrier phase specie i
// - fully explicit
inline tmp<DimensionedField<scalar, volMesh> > Srho
(
const label i
) const;
//- Return total mass source term [kg/m3/s]
inline tmp<fvScalarMatrix> Srho(volScalarField& rho) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "coalCloudListI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,262 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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 "fvMatrices.H"
#include "volFields.H"
#include "DimensionedField.H"
Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
Foam::coalCloudList::UTrans() const
{
tmp<DimensionedField<vector, volMesh> > tfld
(
new DimensionedField<vector, volMesh>
(
IOobject
(
"UTransEff",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedVector("zero", dimMass*dimVelocity, vector::zero)
)
);
DimensionedField<vector, volMesh>& fld = tfld();
forAll(*this, i)
{
fld += operator[](i).UTrans();
}
return tfld;
}
Foam::tmp<Foam::fvVectorMatrix> Foam::coalCloudList::SU
(
volVectorField& U
) const
{
tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
fvVectorMatrix& fvm = tfvm();
forAll(*this, i)
{
fvm += operator[](i).SU(U);
}
return tfvm;
}
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::coalCloudList::hsTrans() const
{
tmp<DimensionedField<scalar, volMesh> > tfld
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
"hsTransEff",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("zero", dimEnergy, 0.0)
)
);
DimensionedField<scalar, volMesh>& fld = tfld();
forAll(*this, i)
{
fld += operator[](i).hsTrans();
}
return tfld;
}
Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Sh
(
volScalarField& hs
) const
{
tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(hs, dimEnergy/dimTime));
fvScalarMatrix& fvm = tfvm();
forAll(*this, i)
{
fvm += operator[](i).Sh(hs);
}
return tfvm;
}
Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::SYi
(
const label ii,
volScalarField& Yi
) const
{
tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(Yi, dimMass/dimTime));
fvScalarMatrix& fvm = tfvm();
forAll(*this, i)
{
fvm += operator[](i).SYi(ii, Yi);
}
return tfvm;
}
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::coalCloudList::rhoTrans() const
{
tmp<DimensionedField<scalar, volMesh> > tfld
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
"rhoTransEff",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("zero", dimMass, 0.0)
)
);
DimensionedField<scalar, volMesh>& fld = tfld();
forAll(*this, i)
{
forAll(operator[](i).rhoTrans(), j)
{
fld += operator[](i).rhoTrans()[j];
}
}
return tfld;
}
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::coalCloudList::Srho() const
{
tmp<DimensionedField<scalar, volMesh> > tfld
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
"rhoTransEff",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0.0)
)
);
DimensionedField<scalar, volMesh>& fld = tfld();
forAll(*this, i)
{
fld += operator[](i).Srho();
}
return tfld;
}
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
Foam::coalCloudList::Srho
(
const label i
) const
{
tmp<DimensionedField<scalar, volMesh> > tfld
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
"rhoTransEff",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("zero", dimDensity/dimTime, 0.0)
)
);
DimensionedField<scalar, volMesh>& fld = tfld();
forAll(*this, j)
{
fld += operator[](j).Srho(i);
}
return tfld;
}
Foam::tmp<Foam::fvScalarMatrix> Foam::coalCloudList::Srho
(
volScalarField& rho
) const
{
tmp<fvScalarMatrix> tfvm(new fvScalarMatrix(rho, dimMass/dimTime));
fvScalarMatrix& fvm = tfvm();
forAll(*this, i)
{
fvm += operator[](i).Srho(rho);
}
return tfvm;
}
// ************************************************************************* //

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -31,6 +31,7 @@ License
#include "NoSurfaceReaction.H"
#include "COxidationDiffusionLimitedRate.H"
#include "COxidationKineticDiffusionLimitedRate.H"
#include "COxidationHurtMitchell.H"
#include "COxidationMurphyShaddix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -43,6 +44,7 @@ License
COxidationKineticDiffusionLimitedRate, \
CloudType \
); \
makeSurfaceReactionModelType(COxidationHurtMitchell, CloudType); \
makeSurfaceReactionModelType(COxidationMurphyShaddix, CloudType);

View File

@ -0,0 +1,206 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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 "COxidationHurtMitchell.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class CloudType>
Foam::COxidationHurtMitchell<CloudType>::COxidationHurtMitchell
(
const dictionary& dict,
CloudType& owner
)
:
SurfaceReactionModel<CloudType>(dict, owner, typeName),
Sb_(readScalar(this->coeffDict().lookup("Sb"))),
CsLocalId_(-1),
ashLocalId_(-1),
O2GlobalId_(owner.composition().globalCarrierId("O2")),
CO2GlobalId_(owner.composition().globalCarrierId("CO2")),
WC_(0.0),
WO2_(0.0),
HcCO2_(0.0),
heatOfReaction_(-1.0)
{
// Determine Cs and ash ids
label idSolid = owner.composition().idSolid();
CsLocalId_ = owner.composition().localId(idSolid, "C");
ashLocalId_ = owner.composition().localId(idSolid, "ash", true);
// Set local copies of thermo properties
WO2_ = owner.thermo().carrier().W(O2GlobalId_);
const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
WC_ = WCO2 - WO2_;
HcCO2_ = owner.thermo().carrier().Hc(CO2GlobalId_);
const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
Info<< " C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
if (this->coeffDict().readIfPresent("heatOfReaction", heatOfReaction_))
{
Info<< " Using user specified heat of reaction: "
<< heatOfReaction_ << " [J/kg]" << endl;
}
}
template<class CloudType>
Foam::COxidationHurtMitchell<CloudType>::COxidationHurtMitchell
(
const COxidationHurtMitchell<CloudType>& srm
)
:
SurfaceReactionModel<CloudType>(srm),
Sb_(srm.Sb_),
CsLocalId_(srm.CsLocalId_),
ashLocalId_(srm.ashLocalId_),
O2GlobalId_(srm.O2GlobalId_),
CO2GlobalId_(srm.CO2GlobalId_),
WC_(srm.WC_),
WO2_(srm.WO2_),
HcCO2_(srm.HcCO2_),
heatOfReaction_(srm.heatOfReaction_)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class CloudType>
Foam::COxidationHurtMitchell<CloudType>::~COxidationHurtMitchell()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CloudType>
Foam::scalar Foam::COxidationHurtMitchell<CloudType>::calculate
(
const scalar dt,
const label cellI,
const scalar d,
const scalar T,
const scalar Tc,
const scalar pc,
const scalar rhoc,
const scalar mass,
const scalarField& YGas,
const scalarField& YLiquid,
const scalarField& YSolid,
const scalarField& YMixture,
const scalar N,
scalarField& dMassGas,
scalarField& dMassLiquid,
scalarField& dMassSolid,
scalarField& dMassSRCarrier
) const
{
const label idGas = CloudType::parcelType::GAS;
const label idSolid = CloudType::parcelType::SLD;
const scalar Ychar = YMixture[idSolid]*YSolid[CsLocalId_];
// Surface combustion until combustible fraction is consumed
if (Ychar < SMALL)
{
return 0.0;
}
const SLGThermo& thermo = this->owner().thermo();
// Local mass fraction of O2 in the carrier phase
const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
// No combustion if no oxygen present
if (YO2 < SMALL)
{
return 0.0;
}
// Conversion from [g/cm^2) to [kg/m^2]
const scalar convSI = 1000.0/10000.0;
// Universal gas constant in [kcal/mol/K]
const scalar RRcal = 1985.877534;
// Dry mass fraction
scalar Ydaf = YMixture[idGas] + YMixture[idSolid];
if (ashLocalId_ != -1)
{
Ydaf -= YMixture[idSolid]*YSolid[ashLocalId_];
}
// Char percentage
const scalar charPrc = Ychar/Ydaf*100.0;
// Particle surface area
const scalar Ap = constant::mathematical::pi*sqr(d);
// Far field partial pressure O2 [Pa]
// Note: Should really use the surface partial pressure
const scalar ppO2 = max(0.0, rhoc*YO2/WO2_*specie::RR*Tc);
// Activation energy [kcal/mol]
const scalar E = -5.94 + 0.355*charPrc;
// Pre-exponential factor [g/(cm^2.s.atm^0.5)]
const scalar lnK1750 = 2.8 - 0.0758*charPrc;
const scalar A = exp(lnK1750 + E/RRcal/1750.0);
// Kinetic rate of char oxidation [g/(cm^2.s.atm^0.5)]
const scalar Rk = A*exp(-E/(RRcal*T));
// Molar reaction rate per unit surface area [kmol/(m^2.s)]
const scalar qCsLim = mass*Ychar/(WC_*Ap*dt);
const scalar qCs = min(convSI*Rk*Foam::sqrt(ppO2/101325.0), qCsLim);
// Calculate the number of molar units reacted [kmol]
const scalar dOmega = qCs*Ap*dt;
// Add to carrier phase mass transfer
dMassSRCarrier[O2GlobalId_] += -dOmega*Sb_*WO2_;
dMassSRCarrier[CO2GlobalId_] += dOmega*(WC_ + Sb_*WO2_);
// Add to particle mass transfer
dMassSolid[CsLocalId_] += dOmega*WC_;
// Return the heat of reaction [J]
// note: carrier sensible enthalpy exchange handled via change in mass
if (heatOfReaction_ < 0)
{
const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
return dOmega*(WC_*HsC - (WC_ + Sb_*WO2_)*HcCO2_);
}
else
{
return dOmega*WC_*heatOfReaction_;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,184 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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
COxidationHurtMitchell
Description
Char oxidation model given by Hurt and Mitchell:
Based on the reference:
Hurt R. and Mitchell R., "Unified high-temperature char combustion
kinetics for a suite of coals of various rank", 24th Symposium in
Combustion, The Combustion Institute, 1992, p 1243-1250
Model specifies the rate of char combustion.
C(s) + Sb*O2 -> CO2
where Sb is the stoichiometry of the reaction
Model validity:
Gas temperature: Tc > 1500 K
Particle sizes: 75 um -> 200 um
Pox > 0.3 atm
\*---------------------------------------------------------------------------*/
#ifndef COxidationHurtMitchell_H
#define COxidationHurtMitchell_H
#include "SurfaceReactionModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward class declarations
template<class CloudType>
class COxidationHurtMitchell;
/*---------------------------------------------------------------------------*\
Class COxidationHurtMitchell Declaration
\*---------------------------------------------------------------------------*/
template<class CloudType>
class COxidationHurtMitchell
:
public SurfaceReactionModel<CloudType>
{
// Private data
// Model constants
//- Stoichiometry of reaction
const scalar Sb_;
// Addressing
//- Cs position in global/local lists
label CsLocalId_;
//- Ash position in global/local lists
label ashLocalId_;
//- O2 position in global list
label O2GlobalId_;
//- CO2 positions in global list
label CO2GlobalId_;
// Local copies of thermo properties
//- Molecular weight of C [kg/kmol]
scalar WC_;
//- Molecular weight of O2 [kg/kmol]
scalar WO2_;
//- Formation enthalpy for CO2 [J/kg]
scalar HcCO2_;
//- Heat of reaction [J/kg] (optional)
scalar heatOfReaction_;
public:
//- Runtime type information
TypeName("COxidationHurtMitchell");
// Constructors
//- Construct from dictionary
COxidationHurtMitchell
(
const dictionary& dict,
CloudType& owner
);
//- Construct copy
COxidationHurtMitchell
(
const COxidationHurtMitchell<CloudType>& srm
);
//- Construct and return a clone
virtual autoPtr<SurfaceReactionModel<CloudType> > clone() const
{
return autoPtr<SurfaceReactionModel<CloudType> >
(
new COxidationHurtMitchell<CloudType>(*this)
);
}
//- Destructor
virtual ~COxidationHurtMitchell();
// Member Functions
//- Update surface reactions
virtual scalar calculate
(
const scalar dt,
const label cellI,
const scalar d,
const scalar T,
const scalar Tc,
const scalar pc,
const scalar rhoc,
const scalar mass,
const scalarField& YGas,
const scalarField& YLiquid,
const scalarField& YSolid,
const scalarField& YMixture,
const scalar N,
scalarField& dMassGas,
scalarField& dMassLiquid,
scalarField& dMassSolid,
scalarField& dMassSRCarrier
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "COxidationHurtMitchell.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -661,7 +661,7 @@ void Foam::KinematicCloud<CloudType>::evolve()
template<class CloudType>
template<class TrackData>
void Foam::KinematicCloud<CloudType>::motion(TrackData& td)
void Foam::KinematicCloud<CloudType>::motion(TrackData& td)
{
td.part() = TrackData::tpLinearTrack;
CloudType::move(td, solution_.trackTime());
@ -670,6 +670,152 @@ void Foam::KinematicCloud<CloudType>::motion(TrackData& td)
}
template<class CloudType>
void Foam::KinematicCloud<CloudType>::patchData
(
const parcelType& p,
const polyPatch& pp,
const scalar trackFraction,
const tetIndices& tetIs,
vector& nw,
vector& Up
) const
{
label patchI = pp.index();
label patchFaceI = pp.whichFace(p.face());
vector n = tetIs.faceTri(mesh_).normal();
n /= mag(n);
vector U = U_.boundaryField()[patchI][patchFaceI];
// Unless the face is rotating, the required normal is n;
nw = n;
if (!mesh_.moving())
{
// Only wall patches may have a non-zero wall velocity from
// the velocity field when the mesh is not moving.
if (isA<wallPolyPatch>(pp))
{
Up = U;
}
else
{
Up = vector::zero;
}
}
else
{
vector U00 = U_.oldTime().boundaryField()[patchI][patchFaceI];
vector n00 = tetIs.oldFaceTri(mesh_).normal();
// Difference in normal over timestep
vector dn = vector::zero;
if (mag(n00) > SMALL)
{
// If the old normal is zero (for example in layer
// addition) then use the current normal, meaning that the
// motion can only be translational, and dn remains zero,
// otherwise, calculate dn:
n00 /= mag(n00);
dn = n - n00;
}
// Total fraction through the timestep of the motion,
// including stepFraction before the current tracking step
// and the current trackFraction
// i.e.
// let s = stepFraction, t = trackFraction
// Motion of x in time:
// |-----------------|---------|---------|
// x00 x0 xi x
//
// where xi is the correct value of x at the required
// tracking instant.
//
// x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
//
// i.e. the motion covered by previous tracking portions
// within this timestep, and
//
// xi = x0 + t*(x - x0)
// = t*x + (1 - t)*x0
// = t*x + (1 - t)*(s*x + (1 - s)*x00)
// = (s + t - s*t)*x + (1 - (s + t - s*t))*x00
//
// let m = (s + t - s*t)
//
// xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
//
// In the same form as before.
scalar m =
p.stepFraction()
+ trackFraction
- (p.stepFraction()*trackFraction);
// When the mesh is moving, the velocity field on wall patches
// will contain the velocity associated with the motion of the
// mesh, in which case it is interpolated in time using m.
// For other patches the face velocity will need to be
// reconstructed from the face centre motion.
const vector& Cf = mesh_.faceCentres()[p.face()];
vector Cf00 = mesh_.faces()[p.face()].centre(mesh_.oldPoints());
if (isA<wallPolyPatch>(pp))
{
Up = U00 + m*(U - U00);
}
else
{
Up = (Cf - Cf00)/mesh_.time().deltaTValue();
}
if (mag(dn) > SMALL)
{
// Rotational motion, nw requires interpolation and a
// rotational velocity around face centre correction to Up
// is required.
nw = n00 + m*dn;
// Cf at tracking instant
vector Cfi = Cf00 + m*(Cf - Cf00);
// Normal vector cross product
vector omega = (n00 ^ n);
scalar magOmega = mag(omega);
// magOmega = sin(angle between unit normals)
// Normalise omega vector by magOmega, then multiply by
// angle/dt to give the correct angular velocity vector.
omega *= Foam::asin(magOmega)/(magOmega*mesh_.time().deltaTValue());
// Project position onto face and calculate this position
// relative to the face centre.
vector facePos =
p.position()
- ((p.position() - Cfi) & nw)*nw
- Cfi;
Up += (omega ^ facePos);
}
// No further action is required if the motion is
// translational only, nw and Up have already been set.
}
}
template<class CloudType>
void Foam::KinematicCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
{

View File

@ -554,6 +554,17 @@ public:
template<class TrackData>
void motion(TrackData& td);
//- Calculate the patch normal and velocity to interact with,
// accounting for patch motion if required.
void patchData
(
const parcelType& p,
const polyPatch& pp,
const scalar trackFraction,
const tetIndices& tetIs,
vector& normal,
vector& Up
) const;
// Mapping

View File

@ -87,8 +87,6 @@ void Foam::ReactingCloud<CloudType>::cloudReset(ReactingCloud<CloudType>& c)
compositionModel_.reset(c.compositionModel_.ptr());
phaseChangeModel_.reset(c.phaseChangeModel_.ptr());
dMassPhaseChange_ = c.dMassPhaseChange_;
}
@ -111,8 +109,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
constProps_(this->particleProperties(), this->solution().active()),
compositionModel_(NULL),
phaseChangeModel_(NULL),
rhoTrans_(thermo.carrier().species().size()),
dMassPhaseChange_(0.0)
rhoTrans_(thermo.carrier().species().size())
{
if (this->solution().active())
{
@ -167,8 +164,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
constProps_(c.constProps_),
compositionModel_(c.compositionModel_->clone()),
phaseChangeModel_(c.phaseChangeModel_->clone()),
rhoTrans_(c.rhoTrans_.size()),
dMassPhaseChange_(c.dMassPhaseChange_)
rhoTrans_(c.rhoTrans_.size())
{
forAll(c.rhoTrans_, i)
{
@ -209,8 +205,7 @@ Foam::ReactingCloud<CloudType>::ReactingCloud
compositionModel_(c.compositionModel_->clone()),
// compositionModel_(NULL),
phaseChangeModel_(NULL),
rhoTrans_(0),
dMassPhaseChange_(0.0)
rhoTrans_(0)
{}

View File

@ -121,12 +121,6 @@ protected:
PtrList<DimensionedField<scalar, volMesh> > rhoTrans_;
// Check
//- Total mass transferred to continuous phase via phase change
scalar dMassPhaseChange_;
// Protected Member Functions
// New parcel helper functions

View File

@ -51,7 +51,69 @@ void Foam::ThermoCloud<CloudType>::setModels()
).ptr()
);
this->subModelProperties().lookup("radiation") >> radiation_;
if (this->solution().coupled())
{
this->subModelProperties().lookup("radiation") >> radiation_;
}
if (radiation_)
{
radAreaP_.reset
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
this->name() + "::radAreaP",
this->db().time().timeName(),
this->db(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
this->mesh(),
dimensionedScalar("zero", dimArea, 0.0)
)
);
radT4_.reset
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
this->name() + "::radT4",
this->db().time().timeName(),
this->db(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
this->mesh(),
dimensionedScalar("zero", pow4(dimTemperature), 0.0)
)
);
radAreaPT4_.reset
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
this->name() + "::radAreaPT4",
this->db().time().timeName(),
this->db(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
this->mesh(),
dimensionedScalar
(
"zero",
sqr(dimLength)*pow4(dimTemperature),
0.0
)
)
);
}
}
@ -98,13 +160,16 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
heatTransferModel_(NULL),
TIntegrator_(NULL),
radiation_(false),
radAreaP_(NULL),
radT4_(NULL),
radAreaPT4_(NULL),
hsTrans_
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
this->name() + "hsTrans",
this->name() + "::hsTrans",
this->db().time().timeName(),
this->db(),
IOobject::READ_IF_PRESENT,
@ -120,7 +185,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
(
IOobject
(
this->name() + "hsCoeff",
this->name() + "::hsCoeff",
this->db().time().timeName(),
this->db(),
IOobject::READ_IF_PRESENT,
@ -165,13 +230,16 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
heatTransferModel_(c.heatTransferModel_->clone()),
TIntegrator_(c.TIntegrator_->clone()),
radiation_(c.radiation_),
radAreaP_(NULL),
radT4_(NULL),
radAreaPT4_(NULL),
hsTrans_
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
this->name() + "hsTrans",
this->name() + "::hsTrans",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
@ -187,7 +255,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
(
IOobject
(
this->name() + "hsCoeff",
this->name() + "::hsCoeff",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
@ -197,7 +265,61 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
c.hsCoeff()
)
)
{}
{
if (radiation_)
{
radAreaP_.reset
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
this->name() + "::radAreaP",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
c.radAreaP()
)
);
radT4_.reset
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
this->name() + "::radT4",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
c.radT4()
)
);
radAreaPT4_.reset
(
new DimensionedField<scalar, volMesh>
(
IOobject
(
this->name() + "::radAreaPT4",
this->db().time().timeName(),
this->db(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
c.radAreaPT4()
)
);
}
}
template<class CloudType>
@ -218,6 +340,9 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
heatTransferModel_(NULL),
TIntegrator_(NULL),
radiation_(false),
radAreaP_(NULL),
radT4_(NULL),
radAreaPT4_(NULL),
hsTrans_(NULL),
hsCoeff_(NULL)
{}
@ -285,6 +410,13 @@ void Foam::ThermoCloud<CloudType>::resetSourceTerms()
CloudType::resetSourceTerms();
hsTrans_->field() = 0.0;
hsCoeff_->field() = 0.0;
if (radiation_)
{
radAreaP_->field() = 0.0;
radT4_->field() = 0.0;
radAreaPT4_->field() = 0.0;
}
}
@ -298,6 +430,13 @@ void Foam::ThermoCloud<CloudType>::relaxSources
this->relax(hsTrans_(), cloudOldTime.hsTrans(), "h");
this->relax(hsCoeff_(), cloudOldTime.hsCoeff(), "h");
if (radiation_)
{
this->relax(radAreaP_(), cloudOldTime.radAreaP(), "radiation");
this->relax(radT4_(), cloudOldTime.radT4(), "radiation");
this->relax(radAreaPT4_(), cloudOldTime.radAreaPT4(), "radiation");
}
}
@ -308,6 +447,13 @@ void Foam::ThermoCloud<CloudType>::scaleSources()
this->scale(hsTrans_(), "h");
this->scale(hsCoeff_(), "h");
if (radiation_)
{
this->scale(radAreaP_(), "radiation");
this->scale(radT4_(), "radiation");
this->scale(radAreaPT4_(), "radiation");
}
}

View File

@ -132,6 +132,15 @@ protected:
//- Include radiation
Switch radiation_;
//- Radiation sum of parcel projected areas
autoPtr<DimensionedField<scalar, volMesh> > radAreaP_;
//- Radiation sum of parcel temperature^4
autoPtr<DimensionedField<scalar, volMesh> > radT4_;
//- Radiation sum of parcel projected areas * temperature^4
autoPtr<DimensionedField<scalar, volMesh> > radAreaPT4_;
// Sources
@ -244,6 +253,26 @@ public:
//- Radiation flag
inline bool radiation() const;
//- Radiation sum of parcel projected areas [m2]
inline DimensionedField<scalar, volMesh>& radAreaP();
//- Radiation sum of parcel projected areas [m2]
inline const DimensionedField<scalar, volMesh>&
radAreaP() const;
//- Radiation sum of parcel temperature^4 [K4]
inline DimensionedField<scalar, volMesh>& radT4();
//- Radiation sum of parcel temperature^4 [K4]
inline const DimensionedField<scalar, volMesh>& radT4() const;
//- Radiation sum of parcel projected area*temperature^4 [m2K4]
inline DimensionedField<scalar, volMesh>& radAreaPT4();
//- Radiation sum of parcel temperature^4 [m2K4]
inline const DimensionedField<scalar, volMesh>&
radAreaPT4() const;
// Sources

View File

@ -89,6 +89,114 @@ inline bool Foam::ThermoCloud<CloudType>::radiation() const
}
template<class CloudType>
inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ThermoCloud<CloudType>::radAreaP()
{
if (!radiation_)
{
FatalErrorIn
(
"inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
"Foam::ThermoCloud<CloudType>::radAreaP()"
) << "Radiation field requested, but radiation model not active"
<< abort(FatalError);
}
return radAreaP_();
}
template<class CloudType>
inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ThermoCloud<CloudType>::radAreaP() const
{
if (!radiation_)
{
FatalErrorIn
(
"inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
"Foam::ThermoCloud<CloudType>::radAreaP()"
) << "Radiation field requested, but radiation model not active"
<< abort(FatalError);
}
return radAreaP_();
}
template<class CloudType>
inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ThermoCloud<CloudType>::radT4()
{
if (!radiation_)
{
FatalErrorIn
(
"inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
"Foam::ThermoCloud<CloudType>::radT4()"
) << "Radiation field requested, but radiation model not active"
<< abort(FatalError);
}
return radT4_();
}
template<class CloudType>
inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ThermoCloud<CloudType>::radT4() const
{
if (!radiation_)
{
FatalErrorIn
(
"inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
"Foam::ThermoCloud<CloudType>::radT4()"
) << "Radiation field requested, but radiation model not active"
<< abort(FatalError);
}
return radT4_();
}
template<class CloudType>
inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ThermoCloud<CloudType>::radAreaPT4()
{
if (!radiation_)
{
FatalErrorIn
(
"inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
"Foam::ThermoCloud<CloudType>::radAreaPT4()"
) << "Radiation field requested, but radiation model not active"
<< abort(FatalError);
}
return radAreaPT4_();
}
template<class CloudType>
inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ThermoCloud<CloudType>::radAreaPT4() const
{
if (!radiation_)
{
FatalErrorIn
(
"inline Foam::DimensionedField<Foam::scalar, Foam::volMesh> "
"Foam::ThermoCloud<CloudType>::radAreaPT4()"
) << "Radiation field requested, but radiation model not active"
<< abort(FatalError);
}
return radAreaPT4_();
}
template<class CloudType>
inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
Foam::ThermoCloud<CloudType>::hsTrans()
@ -182,21 +290,15 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::Ep() const
)
);
// Need to check if coupled as field is created on-the-fly
if (radiation_ && this->solution().coupled())
if (radiation_)
{
scalarField& Ep = tEp().internalField();
const scalar dt = this->db().time().deltaTValue();
const scalarField& V = this->mesh().V();
const scalar epsilon = constProps_.epsilon0();
const scalarField& sumAreaPT4 = radAreaPT4_->field();
forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
{
const parcelType& p = iter();
const label cellI = p.cell();
Ep[cellI] += p.nParticle()*p.areaP()*pow4(p.T());
}
Ep *= epsilon*physicoChemical::sigma.value()/V;
Ep = sumAreaPT4*epsilon*physicoChemical::sigma.value()/V/dt;
}
return tEp;
@ -224,21 +326,15 @@ inline Foam::tmp<Foam::volScalarField> Foam::ThermoCloud<CloudType>::ap() const
)
);
// Need to check if coupled as field is created on-the-fly
if (radiation_ && this->solution().coupled())
if (radiation_)
{
scalarField& ap = tap().internalField();
const scalar dt = this->db().time().deltaTValue();
const scalarField& V = this->mesh().V();
const scalar epsilon = constProps_.epsilon0();
const scalarField& sumAreaP = radAreaP_->field();
forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
{
const parcelType& p = iter();
const label cellI = p.cell();
ap[cellI] += p.nParticle()*p.areaP();
}
ap *= epsilon/V;
ap = sumAreaP*epsilon/V/dt;
}
return tap;
@ -267,23 +363,16 @@ Foam::ThermoCloud<CloudType>::sigmap() const
)
);
// Need to check if coupled as field is created on-the-fly
if (radiation_ && this->solution().coupled())
if (radiation_)
{
scalarField& sigmap = tsigmap().internalField();
const scalar dt = this->db().time().deltaTValue();
const scalarField& V = this->mesh().V();
const scalar epsilon = constProps_.epsilon0();
const scalar f = constProps_.f0();
const scalarField& sumAreaP = radAreaP_->field();
forAllConstIter(typename ThermoCloud<CloudType>, *this, iter)
{
const parcelType& p = iter();
const label cellI = p.cell();
sigmap[cellI] += p.nParticle()*p.areaP();
}
sigmap *= (1.0 - f)*(1.0 - epsilon)/V;
sigmap *= sumAreaP*(1.0 - f)*(1.0 - epsilon)/V/dt;
}
return tsigmap;

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -363,7 +363,7 @@ bool Foam::KinematicParcel<ParcelType>::hitPatch
static_cast<typename TrackData::cloudType::parcelType&>(*this);
// Invoke post-processing model
td.cloud().functions().postPatch(p, patchI, pp.whichFace(p.face()));
td.cloud().functions().postPatch(p, pp, trackFraction, tetIs);
// Invoke surface film model
if (td.cloud().surfaceFilm().transferParcel(p, pp, td.keepParticle))

View File

@ -268,12 +268,15 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
(
td,
dt,
this->age_,
Ts,
d0,
T0,
mass0,
this->mass0_,
YMix[GAS]*YGas_,
YMix[LIQ]*YLiquid_,
YMix[SLD]*YSolid_,
canCombust_,
dMassDV,
Sh,
@ -466,6 +469,16 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// Update sensible enthalpy transfer
td.cloud().hsTrans()[cellI] += np0*dhsTrans;
td.cloud().hsCoeff()[cellI] += np0*Sph;
// Update radiation fields
if (td.cloud().radiation())
{
const scalar ap = this->areaP();
const scalar T4 = pow4(this->T_);
td.cloud().radAreaP()[cellI] += dt*np0*ap;
td.cloud().radT4()[cellI] += dt*np0*T4;
td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
}
}
}
@ -476,12 +489,15 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
(
TrackData& td,
const scalar dt,
const scalar age,
const scalar Ts,
const scalar d,
const scalar T,
const scalar mass,
const scalar mass0,
const scalarField& YGasEff,
const scalarField& YLiquidEff,
const scalarField& YSolidEff,
bool& canCombust,
scalarField& dMassDV,
scalar& Sh,
@ -510,10 +526,13 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
td.cloud().devolatilisation().calculate
(
dt,
age,
mass0,
mass,
T,
YGasEff,
YLiquidEff,
YSolidEff,
canCombust,
dMassDV
);
@ -638,7 +657,8 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
ParcelType(p),
YGas_(p.YGas_),
YLiquid_(p.YLiquid_),
YSolid_(p.YSolid_)
YSolid_(p.YSolid_),
canCombust_(p.canCombust_)
{}
@ -652,7 +672,8 @@ Foam::ReactingMultiphaseParcel<ParcelType>::ReactingMultiphaseParcel
ParcelType(p, mesh),
YGas_(p.YGas_),
YLiquid_(p.YLiquid_),
YSolid_(p.YSolid_)
YSolid_(p.YSolid_),
canCombust_(p.canCombust_)
{}

View File

@ -196,12 +196,15 @@ protected:
(
TrackData& td,
const scalar dt, // timestep
const scalar Ts, // Surface temperature
const scalar age, // age
const scalar Ts, // surface temperature
const scalar d, // diameter
const scalar T, // temperature
const scalar mass, // mass
const scalar mass0, // mass (initial on injection)
const scalarField& YGasEff,// Gas component mass fractions
const scalarField& YGasEff,// gas component mass fractions
const scalarField& YLiquidEff,// liquid component mass fractions
const scalarField& YSolidEff,// solid component mass fractions
bool& canCombust, // 'can combust' flag
scalarField& dMassDV, // mass transfer - local to particle
scalar& Sh, // explicit particle enthalpy source

View File

@ -458,6 +458,16 @@ void Foam::ReactingParcel<ParcelType>::calc
// Update sensible enthalpy transfer
td.cloud().hsTrans()[cellI] += np0*dhsTrans;
td.cloud().hsCoeff()[cellI] += np0*Sph;
// Update radiation fields
if (td.cloud().radiation())
{
const scalar ap = this->areaP();
const scalar T4 = pow4(this->T_);
td.cloud().radAreaP()[cellI] += dt*np0*ap;
td.cloud().radT4()[cellI] += dt*np0*T4;
td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
}
}
}

View File

@ -144,7 +144,7 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
Ts = td.cloud().constProps().TMin();
}
// Assuming thermo props vary linearly with T for small dT
// Assuming thermo props vary linearly with T for small d(T)
const scalar TRatio = Tc_/Ts;
rhos = this->rhoc_*TRatio;
@ -252,6 +252,16 @@ void Foam::ThermoParcel<ParcelType>::calc
// Update sensible enthalpy coefficient
td.cloud().hsCoeff()[cellI] += np0*Sph;
// Update radiation fields
if (td.cloud().radiation())
{
const scalar ap = this->areaP();
const scalar T4 = pow4(this->T_);
td.cloud().radAreaP()[cellI] += dt*np0*ap;
td.cloud().radT4()[cellI] += dt*np0*T4;
td.cloud().radAreaP()[cellI] += dt*np0*ap*T4;
}
}
}

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -107,8 +107,9 @@ template<class CloudType>
void Foam::CloudFunctionObject<CloudType>::postPatch
(
const typename CloudType::parcelType&,
const label,
const label
const polyPatch&,
const scalar,
const tetIndices&
)
{
// do nothing

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -141,8 +141,9 @@ public:
virtual void postPatch
(
const typename CloudType::parcelType& p,
const label patchI,
const label patchFaceI
const polyPatch& pp,
const scalar trackFraction,
const tetIndices& testIs
);
//- Post-face hook

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -146,13 +146,14 @@ template<class CloudType>
void Foam::CloudFunctionObjectList<CloudType>::postPatch
(
const typename CloudType::parcelType& p,
const label patchI,
const label patchFaceI
const polyPatch& pp,
const scalar trackFraction,
const tetIndices& tetIs
)
{
forAll(*this, i)
{
this->operator[](i).postPatch(p, patchI, patchFaceI);
this->operator[](i).postPatch(p, pp, trackFraction, tetIs);
}
}

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -121,8 +121,9 @@ public:
virtual void postPatch
(
const typename CloudType::parcelType& p,
const label patchI,
const label patchFaceI
const polyPatch& pp,
const scalar trackFraction,
const tetIndices& tetIs
);
//- Post-face hook

View File

@ -166,29 +166,32 @@ template<class CloudType>
void Foam::ParticleErosion<CloudType>::postPatch
(
const parcelType& p,
const label patchI,
const label patchFaceI
const polyPatch& pp,
const scalar trackFraction,
const tetIndices& tetIs
)
{
const label patchI = pp.index();
const label localPatchI = applyToPatch(patchI);
if (localPatchI != -1)
{
const fvMesh& mesh = this->owner().mesh();
vector nw;
vector Up;
// patch-normal direction
vector nw = p.currentTetIndices().faceTri(mesh).normal();
this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
// particle direction of travel
const vector& U = p.U();
// particle velocity reletive to patch
const vector& U = p.U() - Up;
// quick reject if particle travelling away from the patch
if ((-nw & U) < 0)
if ((nw & U) < 0)
{
return;
}
nw /= mag(nw);
const scalar magU = mag(U);
const vector Udir = U/magU;
@ -197,6 +200,7 @@ void Foam::ParticleErosion<CloudType>::postPatch
const scalar coeff = p.nParticle()*p.mass()*sqr(magU)/(p_*psi_*K_);
const label patchFaceI = pp.whichFace(p.face());
scalar& Q = QPtr_->boundaryField()[patchI][patchFaceI];
if (tan(alpha) < K_/6.0)
{

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -126,8 +126,9 @@ public:
virtual void postPatch
(
const parcelType& p,
const label patchI,
const label patchFaceI
const polyPatch& pp,
const scalar trackFraction,
const tetIndices& tetIs
);
};

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -218,11 +218,14 @@ template<class CloudType>
void Foam::PatchPostProcessing<CloudType>::postPatch
(
const parcelType& p,
const label patchI,
const label
const polyPatch& pp,
const scalar,
const tetIndices& tetIs
)
{
const label patchI = pp.index();
const label localPatchI = applyToPatch(patchI);
if (localPatchI != -1 && patchData_[localPatchI].size() < maxStoredParcels_)
{
times_[localPatchI].append(this->owner().time().value());

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -127,8 +127,9 @@ public:
virtual void postPatch
(
const parcelType& p,
const label patchI,
const label patchFaceI
const polyPatch& pp,
const scalar trackFraction,
const tetIndices& tetIs
);
};

View File

@ -306,8 +306,9 @@ void Foam::PairCollision<CloudType>::wallInteraction()
this->owner().functions().postPatch
(
p,
patchI,
patchFaceI
mesh.boundaryMesh()[patchI],
1.0,
p.currentTetIndices()
);
}
}

View File

@ -284,6 +284,7 @@ Foam::InjectionModel<CloudType>::InjectionModel(CloudType& owner)
SOI_(0.0),
volumeTotal_(0.0),
massTotal_(0.0),
massFlowRate_(owner.db().time(), "massFlowRate"),
massInjected_(this->template getModelProperty<scalar>("massInjected")),
nInjections_(this->template getModelProperty<label>("nInjections")),
parcelsAddedTotal_
@ -310,6 +311,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
SOI_(0.0),
volumeTotal_(0.0),
massTotal_(0.0),
massFlowRate_(owner.db().time(), "massFlowRate"),
massInjected_(this->template getModelProperty<scalar>("massInjected")),
nInjections_(this->template getModelProperty<scalar>("nInjections")),
parcelsAddedTotal_
@ -335,7 +337,8 @@ Foam::InjectionModel<CloudType>::InjectionModel
}
else
{
this->coeffDict().lookup("massFlowRate") >> massTotal_;
massFlowRate_.reset(this->coeffDict());
massTotal_ = massFlowRate_.value(owner.db().time().value());
}
const word parcelBasisType = this->coeffDict().lookup("parcelBasisType");
@ -384,6 +387,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
SOI_(im.SOI_),
volumeTotal_(im.volumeTotal_),
massTotal_(im.massTotal_),
massFlowRate_(im.massFlowRate_),
massInjected_(im.massInjected_),
nInjections_(im.nInjections_),
parcelsAddedTotal_(im.parcelsAddedTotal_),
@ -458,7 +462,16 @@ Foam::scalar Foam::InjectionModel<CloudType>::volumeToInject
template<class CloudType>
Foam::scalar Foam::InjectionModel<CloudType>::averageParcelMass()
{
label nTotal = parcelsToInject(0.0, timeEnd() - timeStart());
label nTotal = 0.0;
if (this->owner().solution().transient())
{
nTotal = parcelsToInject(0.0, timeEnd() - timeStart());
}
else
{
nTotal = parcelsToInject(0.0, 1.0);
}
return massTotal_/nTotal;
}
@ -594,6 +607,8 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
const polyMesh& mesh = this->owner().mesh();
typename TrackData::cloudType& cloud = td.cloud();
massTotal_ = massFlowRate_.value(mesh.time().value());
// Reset counters
time0_ = 0.0;
label parcelsAdded = 0;

View File

@ -53,6 +53,7 @@ SourceFiles
#include "runTimeSelectionTables.H"
#include "SubModelBase.H"
#include "vector.H"
#include "TimeDataEntry.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -102,6 +103,9 @@ protected:
//- Total mass to inject [kg]
scalar massTotal_;
//- Mass flow rate profile for steady calculations
TimeDataEntry<scalar> massFlowRate_;
//- Total mass injected to date [kg]
scalar massInjected_;

View File

@ -233,7 +233,7 @@ bool Foam::LocalInteraction<CloudType>::correct
vector nw;
vector Up;
this->patchData(p, pp, trackFraction, tetIs, nw, Up);
this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
// Calculate motion relative to patch velocity
U -= Up;

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -179,159 +179,6 @@ bool Foam::PatchInteractionModel<CloudType>::correct
}
template<class CloudType>
void Foam::PatchInteractionModel<CloudType>::patchData
(
typename CloudType::parcelType& p,
const polyPatch& pp,
const scalar trackFraction,
const tetIndices& tetIs,
vector& nw,
vector& Up
) const
{
const fvMesh& mesh = this->owner().mesh();
const volVectorField& Ufield =
mesh.objectRegistry::lookupObject<volVectorField>(UName_);
label patchI = pp.index();
label patchFaceI = pp.whichFace(p.face());
vector n = tetIs.faceTri(mesh).normal();
n /= mag(n);
vector U = Ufield.boundaryField()[patchI][patchFaceI];
// Unless the face is rotating, the required normal is n;
nw = n;
if (!mesh.moving())
{
// Only wall patches may have a non-zero wall velocity from
// the velocity field when the mesh is not moving.
if (isA<wallPolyPatch>(pp))
{
Up = U;
}
else
{
Up = vector::zero;
}
}
else
{
vector U00 = Ufield.oldTime().boundaryField()[patchI][patchFaceI];
vector n00 = tetIs.oldFaceTri(mesh).normal();
// Difference in normal over timestep
vector dn = vector::zero;
if (mag(n00) > SMALL)
{
// If the old normal is zero (for example in layer
// addition) then use the current normal, meaning that the
// motion can only be translational, and dn remains zero,
// otherwise, calculate dn:
n00 /= mag(n00);
dn = n - n00;
}
// Total fraction thought the timestep of the motion,
// including stepFraction before the current tracking step
// and the current trackFraction
// i.e.
// let s = stepFraction, t = trackFraction
// Motion of x in time:
// |-----------------|---------|---------|
// x00 x0 xi x
//
// where xi is the correct value of x at the required
// tracking instant.
//
// x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
//
// i.e. the motion covered by previous tracking portions
// within this timestep, and
//
// xi = x0 + t*(x - x0)
// = t*x + (1 - t)*x0
// = t*x + (1 - t)*(s*x + (1 - s)*x00)
// = (s + t - s*t)*x + (1 - (s + t - s*t))*x00
//
// let m = (s + t - s*t)
//
// xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
//
// In the same form as before.
scalar m =
p.stepFraction()
+ trackFraction
- (p.stepFraction()*trackFraction);
// When the mesh is moving, the velocity field on wall patches
// will contain the velocity associated with the motion of the
// mesh, in which case it is interpolated in time using m.
// For other patches the face velocity will need to be
// reconstructed from the face centre motion.
const vector& Cf = mesh.faceCentres()[p.face()];
vector Cf00 = mesh.faces()[p.face()].centre(mesh.oldPoints());
if (isA<wallPolyPatch>(pp))
{
Up = U00 + m*(U - U00);
}
else
{
Up = (Cf - Cf00)/this->owner().time().deltaTValue();
}
if (mag(dn) > SMALL)
{
// Rotational motion, nw requires interpolation and a
// rotational velocity around face centre correction to Up
// is required.
nw = n00 + m*dn;
// Cf at tracking instant
vector Cfi = Cf00 + m*(Cf - Cf00);
// Normal vector cross product
vector omega = (n00 ^ n);
scalar magOmega = mag(omega);
// magOmega = sin(angle between unit normals)
// Normalise omega vector by magOmega, then multiply by
// angle/dt to give the correct angular velocity vector.
omega *=
Foam::asin(magOmega)
/(magOmega*this->owner().time().deltaTValue());
// Project position onto face and calculate this position
// relative to the face centre.
vector facePos =
p.position()
- ((p.position() - Cfi) & nw)*nw
- Cfi;
Up += (omega ^ facePos);
}
// No further action is required if the motion is
// translational only, nw and Up have already been set.
}
}
template<class CloudType>
void Foam::PatchInteractionModel<CloudType>::info(Ostream& os)
{

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -164,18 +164,6 @@ public:
const tetIndices& tetIs
);
//- Calculate the patch normal and velocity to interact with,
// accounting for patch motion if required.
void patchData
(
typename CloudType::parcelType& p,
const polyPatch& pp,
const scalar trackFraction,
const tetIndices& tetIs,
vector& normal,
vector& Up
) const;
// I-O

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -74,7 +74,7 @@ bool Foam::Rebound<CloudType>::correct
vector nw;
vector Up;
this->patchData(p, pp, trackFraction, tetIs, nw, Up);
this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
// Calculate motion relative to patch velocity
U -= Up;

View File

@ -148,7 +148,7 @@ bool Foam::StandardWallInteraction<CloudType>::correct
vector nw;
vector Up;
this->patchData(p, pp, trackFraction, tetIs, nw, Up);
this->owner().patchData(p, pp, trackFraction, tetIs, nw, Up);
// Calculate motion relative to patch velocity
U -= Up;

View File

@ -156,18 +156,20 @@ Foam::CompositionModel<CloudType>::componentNames(const label phaseI) const
template<class CloudType>
Foam::label Foam::CompositionModel<CloudType>::globalCarrierId
(
const word& cmptName
const word& cmptName,
const bool allowNotFound
) const
{
label id = thermo_.carrierId(cmptName);
if (id < 0)
if (id < 0 && !allowNotFound)
{
FatalErrorIn
(
"Foam::label Foam::CompositionModel<CloudType>::globalCarrierId"
"("
"const word&"
"const word&, "
"const bool"
") const"
) << "Unable to determine global id for requested component "
<< cmptName << ". Available components are " << nl
@ -182,19 +184,21 @@ template<class CloudType>
Foam::label Foam::CompositionModel<CloudType>::globalId
(
const label phaseI,
const word& cmptName
const word& cmptName,
const bool allowNotFound
) const
{
label id = phaseProps_[phaseI].globalId(cmptName);
if (id < 0)
if (id < 0 && !allowNotFound)
{
FatalErrorIn
(
"Foam::label Foam::CompositionModel<CloudType>::globalId"
"("
"const label, "
"const word&"
"const word&, "
"const bool"
") const"
) << "Unable to determine global id for requested component "
<< cmptName << abort(FatalError);
@ -218,19 +222,21 @@ template<class CloudType>
Foam::label Foam::CompositionModel<CloudType>::localId
(
const label phaseI,
const word& cmptName
const word& cmptName,
const bool allowNotFound
) const
{
label id = phaseProps_[phaseI].id(cmptName);
if (id < 0)
if (id < 0 && !allowNotFound)
{
FatalErrorIn
(
"Foam::label Foam::CompositionModel<CloudType>::localId"
"("
"const label, "
"const word&"
"const word&, "
"const bool"
") const"
) << "Unable to determine local id for component " << cmptName
<< abort(FatalError);
@ -244,12 +250,13 @@ template<class CloudType>
Foam::label Foam::CompositionModel<CloudType>::localToGlobalCarrierId
(
const label phaseI,
const label id
const label id,
const bool allowNotFound
) const
{
label gid = phaseProps_[phaseI].globalCarrierIds()[id];
if (gid < 0)
if (gid < 0 && !allowNotFound)
{
FatalErrorIn
(
@ -257,7 +264,8 @@ Foam::label Foam::CompositionModel<CloudType>::localToGlobalCarrierId
"Foam::CompositionModel<CloudType>::localToGlobalCarrierId"
"("
"const label, "
"const label"
"const label, "
"const bool"
") const"
) << "Unable to determine global carrier id for phase "
<< phaseI << " with local id " << id

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -166,22 +166,37 @@ public:
const wordList& componentNames(const label phaseI) const;
//- Return global id of component cmptName in carrier thermo
label globalCarrierId(const word& cmptName) const;
label globalCarrierId
(
const word& cmptName,
const bool allowNotFound = false
) const;
//- Return global id of component cmptName in phase phaseI
label globalId(const label phaseI, const word& cmptName) const;
label globalId
(
const label phaseI,
const word& cmptName,
const bool allowNotFound = false
) const;
//- Return global ids of for phase phaseI
const labelList& globalIds(const label phaseI) const;
//- Return local id of component cmptName in phase phaseI
label localId(const label phaseI, const word& cmptName) const;
label localId
(
const label phaseI,
const word& cmptName,
const bool allowNotFound = false
) const;
//- Return global carrier id of component given local id
label localToGlobalCarrierId
(
const label phaseI,
const label id
const label id,
const bool allowNotFound = false
) const;
//- Return the list of phase phaseI mass fractions

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -102,10 +102,13 @@ template<class CloudType>
void Foam::ConstantRateDevolatilisation<CloudType>::calculate
(
const scalar dt,
const scalar age,
const scalar mass0,
const scalar mass,
const scalar T,
const scalarField& YGasEff,
const scalarField& YLiquidEff,
const scalarField& YSolidEff,
bool& canCombust,
scalarField& dMassDV
) const

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -104,10 +104,13 @@ public:
virtual void calculate
(
const scalar dt,
const scalar age,
const scalar mass0,
const scalar mass,
const scalar T,
const scalarField& YGasEff,
const scalarField& YLiquidEff,
const scalarField& YSolidEff,
bool& canCombust,
scalarField& dMassDV
) const;

View File

@ -78,6 +78,9 @@ void Foam::DevolatilisationModel<CloudType>::calculate
const scalar,
const scalar,
const scalar,
const scalar,
const scalarField&,
const scalarField&,
const scalarField&,
bool&,
scalarField&
@ -91,6 +94,9 @@ void Foam::DevolatilisationModel<CloudType>::calculate
"const scalar, "
"const scalar, "
"const scalar, "
"const scalar, "
"const scalarField&, "
"const scalarField&, "
"const scalarField&, "
"bool&, "
"scalarField&"

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -47,7 +47,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class DevolatilisationModel Declaration
Class DevolatilisationModel Declaration
\*---------------------------------------------------------------------------*/
template<class CloudType>
@ -126,10 +126,13 @@ public:
virtual void calculate
(
const scalar dt,
const scalar age,
const scalar mass0,
const scalar mass,
const scalar T,
const scalarField& YGasEff,
const scalarField& YLiquidEff,
const scalarField& YSolidEff,
bool& canCombust,
scalarField& dMassDV
) const;

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -71,6 +71,9 @@ void Foam::NoDevolatilisation<CloudType>::calculate
const scalar,
const scalar,
const scalar,
const scalar,
const scalarField&,
const scalarField&,
const scalarField&,
bool& canCombust,
scalarField&

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -85,10 +85,13 @@ public:
virtual void calculate
(
const scalar dt,
const scalar age,
const scalar mass0,
const scalar mass,
const scalar T,
const scalarField& YGasEff,
const scalarField& YLiquidEff,
const scalarField& YSolidEff,
bool& canCombust,
scalarField& dMassDV
) const;

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -105,10 +105,13 @@ template<class CloudType>
void Foam::SingleKineticRateDevolatilisation<CloudType>::calculate
(
const scalar dt,
const scalar age,
const scalar mass0,
const scalar mass,
const scalar T,
const scalarField& YGasEff,
const scalarField& YLiquidEff,
const scalarField& YSolidEff,
bool& canCombust,
scalarField& dMassDV
) const

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-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -229,10 +229,13 @@ public:
virtual void calculate
(
const scalar dt,
const scalar age,
const scalar mass0,
const scalar mass,
const scalar T,
const scalarField& YGasEff,
const scalarField& YLiquidEff,
const scalarField& YSolidEff,
bool& canCombust,
scalarField& dMassDV
) const;

View File

@ -25,7 +25,6 @@ License
#include "AMIInterpolation.H"
#include "meshTools.H"
#include "mergePoints.H"
#include "mapDistribute.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -204,455 +203,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::resetTree
}
template<class SourcePatch, class TargetPatch>
Foam::label Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcDistribution
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch
) const
{
label procI = 0;
if (Pstream::parRun())
{
List<label> facesPresentOnProc(Pstream::nProcs(), 0);
if ((srcPatch.size() > 0) || (tgtPatch.size() > 0))
{
facesPresentOnProc[Pstream::myProcNo()] = 1;
}
else
{
facesPresentOnProc[Pstream::myProcNo()] = 0;
}
Pstream::gatherList(facesPresentOnProc);
Pstream::scatterList(facesPresentOnProc);
label nHaveFaces = sum(facesPresentOnProc);
if (nHaveFaces > 1)
{
procI = -1;
if (debug)
{
Info<< "AMIInterpolation::calcDistribution: "
<< "AMI split across multiple processors" << endl;
}
}
else if (nHaveFaces == 1)
{
procI = findIndex(facesPresentOnProc, 1);
if (debug)
{
Info<< "AMIInterpolation::calcDistribution: "
<< "AMI local to processor" << procI << endl;
}
}
}
// Either not parallel or no faces on any processor
return procI;
}
template<class SourcePatch, class TargetPatch>
Foam::label
Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcOverlappingProcs
(
const List<treeBoundBoxList>& procBb,
const treeBoundBox& bb,
boolList& overlaps
) const
{
overlaps.setSize(procBb.size());
overlaps = false;
label nOverlaps = 0;
forAll(procBb, procI)
{
const List<treeBoundBox>& bbs = procBb[procI];
forAll(bbs, bbI)
{
if (bbs[bbI].overlaps(bb))
{
overlaps[procI] = true;
nOverlaps++;
break;
}
}
}
return nOverlaps;
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::distributePatches
(
const mapDistribute& map,
const TargetPatch& pp,
const globalIndex& gi,
List<faceList>& faces,
List<pointField>& points,
List<labelList>& faceIDs
) const
{
PstreamBuffers pBufs(Pstream::nonBlocking);
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
const labelList& sendElems = map.subMap()[domain];
if (domain != Pstream::myProcNo() && sendElems.size())
{
labelList globalElems(sendElems.size());
forAll(sendElems, i)
{
globalElems[i] = gi.toGlobal(sendElems[i]);
}
faceList subFaces(UIndirectList<face>(pp, sendElems));
primitivePatch subPatch
(
SubList<face>(subFaces, subFaces.size()),
pp.points()
);
if (debug & 2)
{
Pout<< "distributePatches: to processor " << domain
<< " sending faces " << subPatch.faceCentres() << endl;
}
UOPstream toDomain(domain, pBufs);
toDomain
<< subPatch.localFaces() << subPatch.localPoints()
<< globalElems;
}
}
// Start receiving
pBufs.finishedSends();
faces.setSize(Pstream::nProcs());
points.setSize(Pstream::nProcs());
faceIDs.setSize(Pstream::nProcs());
{
// Set up 'send' to myself
const labelList& sendElems = map.subMap()[Pstream::myProcNo()];
faceList subFaces(UIndirectList<face>(pp, sendElems));
primitivePatch subPatch
(
SubList<face>(subFaces, subFaces.size()),
pp.points()
);
// Receive
if (debug & 2)
{
Pout<< "distributePatches: to processor " << Pstream::myProcNo()
<< " sending faces " << subPatch.faceCentres() << endl;
}
faces[Pstream::myProcNo()] = subPatch.localFaces();
points[Pstream::myProcNo()] = subPatch.localPoints();
faceIDs[Pstream::myProcNo()].setSize(sendElems.size());
forAll(sendElems, i)
{
faceIDs[Pstream::myProcNo()][i] = gi.toGlobal(sendElems[i]);
}
}
// Consume
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
const labelList& recvElems = map.constructMap()[domain];
if (domain != Pstream::myProcNo() && recvElems.size())
{
UIPstream str(domain, pBufs);
str >> faces[domain]
>> points[domain]
>> faceIDs[domain];
}
}
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::
distributeAndMergePatches
(
const mapDistribute& map,
const TargetPatch& tgtPatch,
const globalIndex& gi,
faceList& tgtFaces,
pointField& tgtPoints,
labelList& tgtFaceIDs
) const
{
// Exchange per-processor data
List<faceList> allFaces;
List<pointField> allPoints;
List<labelList> allTgtFaceIDs;
distributePatches(map, tgtPatch, gi, allFaces, allPoints, allTgtFaceIDs);
// Renumber and flatten
label nFaces = 0;
label nPoints = 0;
forAll(allFaces, procI)
{
nFaces += allFaces[procI].size();
nPoints += allPoints[procI].size();
}
tgtFaces.setSize(nFaces);
tgtPoints.setSize(nPoints);
tgtFaceIDs.setSize(nFaces);
nFaces = 0;
nPoints = 0;
// My own data first
{
const labelList& faceIDs = allTgtFaceIDs[Pstream::myProcNo()];
SubList<label>(tgtFaceIDs, faceIDs.size()).assign(faceIDs);
const faceList& fcs = allFaces[Pstream::myProcNo()];
forAll(fcs, i)
{
const face& f = fcs[i];
face& newF = tgtFaces[nFaces++];
newF.setSize(f.size());
forAll(f, fp)
{
newF[fp] = f[fp] + nPoints;
}
}
const pointField& pts = allPoints[Pstream::myProcNo()];
forAll(pts, i)
{
tgtPoints[nPoints++] = pts[i];
}
}
// Other proc data follows
forAll(allFaces, procI)
{
if (procI != Pstream::myProcNo())
{
const labelList& faceIDs = allTgtFaceIDs[procI];
SubList<label>(tgtFaceIDs, faceIDs.size(), nFaces).assign(faceIDs);
const faceList& fcs = allFaces[procI];
forAll(fcs, i)
{
const face& f = fcs[i];
face& newF = tgtFaces[nFaces++];
newF.setSize(f.size());
forAll(f, fp)
{
newF[fp] = f[fp] + nPoints;
}
}
const pointField& pts = allPoints[procI];
forAll(pts, i)
{
tgtPoints[nPoints++] = pts[i];
}
}
}
// Merge
labelList oldToNew;
pointField newTgtPoints;
bool hasMerged = mergePoints
(
tgtPoints,
SMALL,
false,
oldToNew,
newTgtPoints
);
if (hasMerged)
{
if (debug)
{
Pout<< "Merged from " << tgtPoints.size()
<< " down to " << newTgtPoints.size() << " points" << endl;
}
tgtPoints.transfer(newTgtPoints);
forAll(tgtFaces, i)
{
inplaceRenumber(oldToNew, tgtFaces[i]);
}
}
}
template<class SourcePatch, class TargetPatch>
Foam::autoPtr<Foam::mapDistribute>
Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcProcMap
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch
) const
{
// Get decomposition of patch
List<treeBoundBoxList> procBb(Pstream::nProcs());
if (srcPatch.size())
{
procBb[Pstream::myProcNo()] = treeBoundBoxList
(
1, // For now single bounding box per proc
treeBoundBox
(
srcPatch.points(),
srcPatch.meshPoints()
)
);
}
else
{
procBb[Pstream::myProcNo()] = treeBoundBoxList();
}
// slightly increase size of bounding boxes to allow for cases where
// bounding boxes are perfectly alligned
forAll(procBb[Pstream::myProcNo()], bbI)
{
treeBoundBox& bb = procBb[Pstream::myProcNo()][bbI];
bb.inflate(0.01);
}
Pstream::gatherList(procBb);
Pstream::scatterList(procBb);
if (debug)
{
Info<< "Determining extent of srcPatch per processor:" << nl
<< "\tproc\tbb" << endl;
forAll(procBb, procI)
{
Info<< '\t' << procI << '\t' << procBb[procI] << endl;
}
}
// Determine which faces of tgtPatch overlaps srcPatch per proc
const faceList& faces = tgtPatch.localFaces();
const pointField& points = tgtPatch.localPoints();
labelListList sendMap;
{
// Per processor indices into all segments to send
List<DynamicList<label> > dynSendMap(Pstream::nProcs());
// Work array - whether processor bb overlaps the face bounds
boolList procBbOverlaps(Pstream::nProcs());
forAll(faces, faceI)
{
if (faces[faceI].size())
{
treeBoundBox faceBb(points, faces[faceI]);
// Find the processor this face overlaps
calcOverlappingProcs(procBb, faceBb, procBbOverlaps);
forAll(procBbOverlaps, procI)
{
if (procBbOverlaps[procI])
{
dynSendMap[procI].append(faceI);
}
}
}
}
// Convert dynamicList to labelList
sendMap.setSize(Pstream::nProcs());
forAll(sendMap, procI)
{
sendMap[procI].transfer(dynSendMap[procI]);
}
}
// Debug printing
if (debug)
{
Pout<< "Of my " << faces.size() << " I need to send to:" << nl
<< "\tproc\tfaces" << endl;
forAll(sendMap, procI)
{
Pout<< '\t' << procI << '\t' << sendMap[procI].size() << endl;
}
}
// Send over how many faces I need to receive
labelListList sendSizes(Pstream::nProcs());
sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
forAll(sendMap, procI)
{
sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
}
Pstream::gatherList(sendSizes);
Pstream::scatterList(sendSizes);
// Determine order of receiving
labelListList constructMap(Pstream::nProcs());
// My local segment first
constructMap[Pstream::myProcNo()] = identity
(
sendMap[Pstream::myProcNo()].size()
);
label segmentI = constructMap[Pstream::myProcNo()].size();
forAll(constructMap, procI)
{
if (procI != Pstream::myProcNo())
{
// What I need to receive is what other processor is sending to me
label nRecv = sendSizes[procI][Pstream::myProcNo()];
constructMap[procI].setSize(nRecv);
for (label i = 0; i < nRecv; i++)
{
constructMap[procI][i] = segmentI++;
}
}
}
autoPtr<mapDistribute> mapPtr
(
new mapDistribute
(
segmentI, // size after construction
sendMap.xfer(),
constructMap.xfer()
)
);
return mapPtr;
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::projectPointsToSurface
(

View File

@ -554,6 +554,7 @@ public:
#ifdef NoRepository
# include "AMIInterpolation.C"
# include "AMIInterpolationParallelOps.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -0,0 +1,481 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "AMIInterpolation.H"
#include "mergePoints.H"
#include "mapDistribute.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::label Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcDistribution
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch
) const
{
label procI = 0;
if (Pstream::parRun())
{
List<label> facesPresentOnProc(Pstream::nProcs(), 0);
if ((srcPatch.size() > 0) || (tgtPatch.size() > 0))
{
facesPresentOnProc[Pstream::myProcNo()] = 1;
}
else
{
facesPresentOnProc[Pstream::myProcNo()] = 0;
}
Pstream::gatherList(facesPresentOnProc);
Pstream::scatterList(facesPresentOnProc);
label nHaveFaces = sum(facesPresentOnProc);
if (nHaveFaces > 1)
{
procI = -1;
if (debug)
{
Info<< "AMIInterpolation::calcDistribution: "
<< "AMI split across multiple processors" << endl;
}
}
else if (nHaveFaces == 1)
{
procI = findIndex(facesPresentOnProc, 1);
if (debug)
{
Info<< "AMIInterpolation::calcDistribution: "
<< "AMI local to processor" << procI << endl;
}
}
}
// Either not parallel or no faces on any processor
return procI;
}
template<class SourcePatch, class TargetPatch>
Foam::label
Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcOverlappingProcs
(
const List<treeBoundBoxList>& procBb,
const treeBoundBox& bb,
boolList& overlaps
) const
{
overlaps.setSize(procBb.size());
overlaps = false;
label nOverlaps = 0;
forAll(procBb, procI)
{
const List<treeBoundBox>& bbs = procBb[procI];
forAll(bbs, bbI)
{
if (bbs[bbI].overlaps(bb))
{
overlaps[procI] = true;
nOverlaps++;
break;
}
}
}
return nOverlaps;
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::distributePatches
(
const mapDistribute& map,
const TargetPatch& pp,
const globalIndex& gi,
List<faceList>& faces,
List<pointField>& points,
List<labelList>& faceIDs
) const
{
PstreamBuffers pBufs(Pstream::nonBlocking);
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
const labelList& sendElems = map.subMap()[domain];
if (domain != Pstream::myProcNo() && sendElems.size())
{
labelList globalElems(sendElems.size());
forAll(sendElems, i)
{
globalElems[i] = gi.toGlobal(sendElems[i]);
}
faceList subFaces(UIndirectList<face>(pp, sendElems));
primitivePatch subPatch
(
SubList<face>(subFaces, subFaces.size()),
pp.points()
);
if (debug & 2)
{
Pout<< "distributePatches: to processor " << domain
<< " sending faces " << subPatch.faceCentres() << endl;
}
UOPstream toDomain(domain, pBufs);
toDomain
<< subPatch.localFaces() << subPatch.localPoints()
<< globalElems;
}
}
// Start receiving
pBufs.finishedSends();
faces.setSize(Pstream::nProcs());
points.setSize(Pstream::nProcs());
faceIDs.setSize(Pstream::nProcs());
{
// Set up 'send' to myself
const labelList& sendElems = map.subMap()[Pstream::myProcNo()];
faceList subFaces(UIndirectList<face>(pp, sendElems));
primitivePatch subPatch
(
SubList<face>(subFaces, subFaces.size()),
pp.points()
);
// Receive
if (debug & 2)
{
Pout<< "distributePatches: to processor " << Pstream::myProcNo()
<< " sending faces " << subPatch.faceCentres() << endl;
}
faces[Pstream::myProcNo()] = subPatch.localFaces();
points[Pstream::myProcNo()] = subPatch.localPoints();
faceIDs[Pstream::myProcNo()].setSize(sendElems.size());
forAll(sendElems, i)
{
faceIDs[Pstream::myProcNo()][i] = gi.toGlobal(sendElems[i]);
}
}
// Consume
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
const labelList& recvElems = map.constructMap()[domain];
if (domain != Pstream::myProcNo() && recvElems.size())
{
UIPstream str(domain, pBufs);
str >> faces[domain]
>> points[domain]
>> faceIDs[domain];
}
}
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::
distributeAndMergePatches
(
const mapDistribute& map,
const TargetPatch& tgtPatch,
const globalIndex& gi,
faceList& tgtFaces,
pointField& tgtPoints,
labelList& tgtFaceIDs
) const
{
// Exchange per-processor data
List<faceList> allFaces;
List<pointField> allPoints;
List<labelList> allTgtFaceIDs;
distributePatches(map, tgtPatch, gi, allFaces, allPoints, allTgtFaceIDs);
// Renumber and flatten
label nFaces = 0;
label nPoints = 0;
forAll(allFaces, procI)
{
nFaces += allFaces[procI].size();
nPoints += allPoints[procI].size();
}
tgtFaces.setSize(nFaces);
tgtPoints.setSize(nPoints);
tgtFaceIDs.setSize(nFaces);
nFaces = 0;
nPoints = 0;
// My own data first
{
const labelList& faceIDs = allTgtFaceIDs[Pstream::myProcNo()];
SubList<label>(tgtFaceIDs, faceIDs.size()).assign(faceIDs);
const faceList& fcs = allFaces[Pstream::myProcNo()];
forAll(fcs, i)
{
const face& f = fcs[i];
face& newF = tgtFaces[nFaces++];
newF.setSize(f.size());
forAll(f, fp)
{
newF[fp] = f[fp] + nPoints;
}
}
const pointField& pts = allPoints[Pstream::myProcNo()];
forAll(pts, i)
{
tgtPoints[nPoints++] = pts[i];
}
}
// Other proc data follows
forAll(allFaces, procI)
{
if (procI != Pstream::myProcNo())
{
const labelList& faceIDs = allTgtFaceIDs[procI];
SubList<label>(tgtFaceIDs, faceIDs.size(), nFaces).assign(faceIDs);
const faceList& fcs = allFaces[procI];
forAll(fcs, i)
{
const face& f = fcs[i];
face& newF = tgtFaces[nFaces++];
newF.setSize(f.size());
forAll(f, fp)
{
newF[fp] = f[fp] + nPoints;
}
}
const pointField& pts = allPoints[procI];
forAll(pts, i)
{
tgtPoints[nPoints++] = pts[i];
}
}
}
// Merge
labelList oldToNew;
pointField newTgtPoints;
bool hasMerged = mergePoints
(
tgtPoints,
SMALL,
false,
oldToNew,
newTgtPoints
);
if (hasMerged)
{
if (debug)
{
Pout<< "Merged from " << tgtPoints.size()
<< " down to " << newTgtPoints.size() << " points" << endl;
}
tgtPoints.transfer(newTgtPoints);
forAll(tgtFaces, i)
{
inplaceRenumber(oldToNew, tgtFaces[i]);
}
}
}
template<class SourcePatch, class TargetPatch>
Foam::autoPtr<Foam::mapDistribute>
Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcProcMap
(
const SourcePatch& srcPatch,
const TargetPatch& tgtPatch
) const
{
// Get decomposition of patch
List<treeBoundBoxList> procBb(Pstream::nProcs());
if (srcPatch.size())
{
procBb[Pstream::myProcNo()] = treeBoundBoxList
(
1, // For now single bounding box per proc
treeBoundBox
(
srcPatch.points(),
srcPatch.meshPoints()
)
);
}
else
{
procBb[Pstream::myProcNo()] = treeBoundBoxList();
}
// slightly increase size of bounding boxes to allow for cases where
// bounding boxes are perfectly alligned
forAll(procBb[Pstream::myProcNo()], bbI)
{
treeBoundBox& bb = procBb[Pstream::myProcNo()][bbI];
bb.inflate(0.01);
}
Pstream::gatherList(procBb);
Pstream::scatterList(procBb);
if (debug)
{
Info<< "Determining extent of srcPatch per processor:" << nl
<< "\tproc\tbb" << endl;
forAll(procBb, procI)
{
Info<< '\t' << procI << '\t' << procBb[procI] << endl;
}
}
// Determine which faces of tgtPatch overlaps srcPatch per proc
const faceList& faces = tgtPatch.localFaces();
const pointField& points = tgtPatch.localPoints();
labelListList sendMap;
{
// Per processor indices into all segments to send
List<DynamicList<label> > dynSendMap(Pstream::nProcs());
// Work array - whether processor bb overlaps the face bounds
boolList procBbOverlaps(Pstream::nProcs());
forAll(faces, faceI)
{
if (faces[faceI].size())
{
treeBoundBox faceBb(points, faces[faceI]);
// Find the processor this face overlaps
calcOverlappingProcs(procBb, faceBb, procBbOverlaps);
forAll(procBbOverlaps, procI)
{
if (procBbOverlaps[procI])
{
dynSendMap[procI].append(faceI);
}
}
}
}
// Convert dynamicList to labelList
sendMap.setSize(Pstream::nProcs());
forAll(sendMap, procI)
{
sendMap[procI].transfer(dynSendMap[procI]);
}
}
// Debug printing
if (debug)
{
Pout<< "Of my " << faces.size() << " I need to send to:" << nl
<< "\tproc\tfaces" << endl;
forAll(sendMap, procI)
{
Pout<< '\t' << procI << '\t' << sendMap[procI].size() << endl;
}
}
// Send over how many faces I need to receive
labelListList sendSizes(Pstream::nProcs());
sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
forAll(sendMap, procI)
{
sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
}
Pstream::gatherList(sendSizes);
Pstream::scatterList(sendSizes);
// Determine order of receiving
labelListList constructMap(Pstream::nProcs());
// My local segment first
constructMap[Pstream::myProcNo()] = identity
(
sendMap[Pstream::myProcNo()].size()
);
label segmentI = constructMap[Pstream::myProcNo()].size();
forAll(constructMap, procI)
{
if (procI != Pstream::myProcNo())
{
// What I need to receive is what other processor is sending to me
label nRecv = sendSizes[procI][Pstream::myProcNo()];
constructMap[procI].setSize(nRecv);
for (label i = 0; i < nRecv; i++)
{
constructMap[procI][i] = segmentI++;
}
}
}
autoPtr<mapDistribute> mapPtr
(
new mapDistribute
(
segmentI, // size after construction
sendMap.xfer(),
constructMap.xfer()
)
);
return mapPtr;
}
// ************************************************************************* //

View File

@ -4,9 +4,11 @@ cd ${0%/*} || exit 1 # run from this directory
# Source tutorial clean functions
. $WM_PROJECT_DIR/bin/tools/CleanFunctions
# remove surface
\rm -f constant/triSurface/propellerTip.obj.gz
cleanCase
\rm -rf 0
\rm -rf constant/extendedFeatureEdgeMesh/
\rm -f constant/triSurface/*.eMesh*